NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


DISSERTATION 


THEORY  OF  MULTIRATE  STATISTICAL 
SIGNAL  PROCESSING  AND 
APPLICATIONS 

by 

Ryan  J.  Kuchler 
September  2005 

Dissertation  Supervisor:  Charles  W.  Therrien 


Approved  for  public  release;  distribution  is  unlimited. 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instruction, 
searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send 
comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing  this  burden, 
to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204, 

Arlington,  Va  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188)  Washington  DC  20503. 

1.  AGENCY  USE  ONLY  (Leave  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 


September  2005  Doctoral  Dissertation 


4.  TITLE  AND  SUBTITLE  Theory  of  Multirate  Statistical 

Signal  Processing  and  Applications 

5.  FUNDING  NUMBERS 

6.  AUTHORS  Kuchler,  Ryan  J. 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Naval  Postgraduate  School 

Monterey  CA  93943-5000 

8.  PERFORMING 

ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 

11.  SUPPLEMENTARY  NOTES  The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect 
the  official  policy  or  position  of  the  Department  of  Defense  or  the  U.S.  Government. 

12a.  DISTRIBUTION/ AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 

12b.  DISTRIBUTION  CODE 

13.  A B ST R A C T ( maximum  200  words) 

This  dissertation  develops  basic  theory  and  applications  of  statistical  multirate  signal  processing.  Specific 
tools  and  terminology  for  describing  multirate  systems  in  the  time  and  frequency  domains  are  presented.  An 
optimal  multirate  estimator  is  derived  in  both  a  direct  form  and  recursive  form.  The  recursive  form  of  the  optimal 
estimator  allows  calculation  of  the  relative  change  in  performance  when  input  signals  are  added  or  removed 
from  the  multirate  system.  The  optimal  multirate  filtering  problem  also  is  specialized  to  the  case  of  optimal 
multirate  linear  prediction.  An  efficient  method  for  calculating  the  multirate  linear  prediction  coefficients  and 
error  variances  is  developed  through  the  use  of  the  multichannel  Levinson  recursion  and  generalized  triangular 
UL  factorization.  Finally,  a  multirate  sequential  classifier  is  derived  and  applied  to  the  problem  of  target 
classification.  It  is  shown  that  classifier  parameters  needed  for  implementing  the  multirate  sequential  classifier 
are  the  same  as  those  for  multirate  linear  prediction.  The  methods  presented  in  this  dissertation  are  useful  for 
multisensor  fusion,  particularly  when  the  sensors  are  operating  at  different  rates. 

14.  SUBJECT  TERMS  Multirate  Statistical  Signal  Processing,  Linear  Estimation,  15.  NUMBER  OF 

Linear  Prediction,  Classification  PAGES  195 

16.  PRICE  CODE 

17.  SECURITY  CLASSIFI-  18.  SECURITY  CLASSIFI-  19.  SECURITY  CLASSIFI-  20.  LIMITATION 

CATION  OF  REPORT  CATION  OF  THIS  PAGE  CATION  OF  ABSTRACT  OF  ABSTRACT 

Unclassified  Unclassified  Unclassified  UL 


NSN  7540-01-280-5500  Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Std.  239-18  298-102 


1 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


n 


Approved  for  public  release;  distribution  is  unlimited 


THEORY  OF  MULTIRATE  STATISTICAL  SIGNAL 

PROCESSING 
AND  APPLICATIONS 

Ryan  J.  Kuchler 

Lieutenant  Commander,  United  States  Navy 
B.S.,  United  States  Naval  Academy,  1993 
M.S.,  Naval  Postgraduate  School,  2002 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 

DOCTOR  OF  PHILOSOPHY  IN  ELECTRICAL  ENGINEERING 

from  the 

NAVAL  POSTGRADUATE  SCHOOL 
September  2005 


Author: 


Ryan  J.  Kuchler 


Approved  by:  Charles  W.  Therrien,  Professor  of  Electrical  Engineering 

Dissertation  Supervisor  and  Committee  Chair 


Roberto  Cristi 
Professor  of  Electrical 
Engineering 

Carlos  F.  Borges 
Associate  Professor  of 
Mathematics 


Murali  Tummala 
Professor  of  Electrical 
Engineering 

Anthony  J.  Healey 
Professor  of  Mechanical 
Engineering 


Approved  by:  Jeffrey  B.  Knorr,  Chair,  Department  of  Electrical  and 

Computer  Engineering 


Approved  by: 


Knox  T.  Millsaps,  Associate  Provost  for  Academic  Affairs 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


IV 


ABSTRACT 


This  dissertation  develops  basic  theory  and  applications  of  statistical  multirate 
signal  processing.  Specific  tools  and  terminology  for  describing  multirate  systems  in 
the  time  and  frequency  domains  are  presented.  An  optimal  multirate  estimator  is 
derived  in  both  a  direct  form  and  recursive  form.  The  recursive  form  of  the  optimal 
estimator  allows  calculation  of  the  relative  change  in  performance  when  input  signals 
are  added  or  removed  from  the  multirate  system.  The  optimal  multirate  filtering 
problem  also  is  specialized  to  the  case  of  optimal  multirate  linear  prediction.  An 
efficient  method  for  calculating  the  multirate  linear  prediction  coefficients  and  error 
variances  is  developed  through  the  use  of  the  multichannel  Levinson  recursion  and 
generalized  triangular  UL  factorization.  Finally,  a  multirate  sequential  classifier  is 
derived  and  applied  to  the  problem  of  target  classification.  It  is  shown  that  classifier 
parameters  needed  for  implementing  the  multirate  sequential  classifier  are  the  same  as 
those  for  multirate  linear  prediction.  The  methods  presented  in  this  dissertation  are 
useful  for  multisensor  fusion,  particularly  when  the  sensors  are  operating  at  different 
rates. 
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EXECUTIVE  SUMMARY 


Sensor  fusion  is  an  important  component  of  target  detection  and  classification. 
The  ability  to  exploit  correlated  information  among  multiple  sensors  with  different 
operating  characteristics  can  lead  to  improved  detection  and  classification  perfor¬ 
mance.  This  dissertation  addresses  multirate  sequential  classification  through  the 
development  of  a  common  multirate  framework  and  the  multirate  linear  prediction 
equations. 

In  order  to  lay  the  foundation  for  representing  multirate  signals,  key  terms 
applicable  to  multirate  systems  are  defined.  These  are  the  fundamental  rate,  system 
period,  system  phase  and  decimation  factors.  These  terms  are  used  to  explicitly 
describe  the  periodic  nature  of  the  multirate  system  and  allow  for  the  appropriate 
selection  of  values  for  the  periodic  components  in  the  multirate  system.  In  addition, 
fundamental  building  blocks  of  the  multirate  system  are  presented,  along  with  the 
statistical  characterizations  of  multirate  signals  as  they  pass  through  the  multirate 
building  blocks.  In  particular,  a  linear  periodically  time-varying  filter  is  presented 
that  permits  the  input  and  output  signals  to  be  sampled  at  different  rates. 

The  first  multirate  optimal  filter  considered  in  this  dissertation  is  the  optimal 
linear  estimator.  The  explicit  direct  form  that  estimates  multiple  input  signals  at 
possibly  different  sampling  rates  is  first  presented  followed  by  a  recursive  innovations 
form.  This  innovations  form  of  the  optimal  estimator  separates  the  direct  form  op¬ 
timal  filters  into  modified  optimal  filters  and  cross  filters.  These  cross  filters  remove 
any  information  from  one  signal  that  is  contained  in  the  other  signals,  in  an  ordered 
fashion.  In  essence,  a  new  set  of  input  signals  that  are  mutually  orthogonal  are  de¬ 
rived  and  used  as  the  inputs  to  the  modified  optimal  filters.  Using  the  innovations 
form  of  the  optimal  filter  allows  one  to  calculate  the  relative  change  in  performance 
in  the  optimal  estimator  when  input  signals  are  added  or  removed  from  the  system. 


xvii 


The  optimal  filtering  problem  is  specialized  to  the  case  of  optimal  multirate 
linear  prediction.  The  multirate  Normal  equations  are  derived  for  a  multirate  system 
with  multiple  input  and  output  signals  which  are  observed  at  different  sampling  rates. 
For  a  multirate  system  with  a  system  period  K,  there  are  up  to  K  distinct  sets  of 
prediction  coefficients  and  error  covariance  matrices  that  apply  in  a  periodic  fashion. 
An  efficient  method  for  calculating  the  multirate  linear  prediction  coefficients  and 
error  variances  is  developed  through  the  use  of  the  multichannel  Levinson  recursion 
and  generalized  triangular  UL  factorization. 

Finally,  a  multirate  sequential  classifier  is  derived  starting  from  the  basic  the¬ 
ory  of  sequential  hypothesis  testing,  ft  is  shown  that  classifier  parameters  needed  for 
implementing  the  multirate  sequential  classifier  are  the  same  as  those  for  multirate 
linear  prediction.  A  multirate  sequential  classifier  is  then  implemented  and  tested  us¬ 
ing  audio  hies  of  a  propellor  plane  and  three  A-10  jet  aircraft.  The  experiments  tested 
the  classifier  performance  in  selecting  between  the  propellor  plane  and  jet  aircraft  as 
certain  system  parameters  are  changed.  These  parameters  are  the  signal-to-noise 
ratios  of  the  observed  signal  and  the  length  of  the  training  data.  In  addition,  the 
performance  of  the  multirate  classifier  is  compared  to  that  of  the  single-channel  and 
multichannel  classifiers  using  similar  data. 
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I.  INTRODUCTION 


Within  the  field  of  signal  processing,  a  growing  interest  in  the  area  of  multirate 
signal  processing  has  developed  over  the  past  three  decades.  This  area  focuses  on 
the  processing  of  systems  that  contain  multiple  signals  that  can  occur  at  different 
sampling  rates.  Since  this  area  deals  with  systems  that  have  components  operating  at 
different  sampling  rates,  processing  techniques  and  descriptions  are  needed  to  account 
for  any  disparity  between  these  sampling  rates.  Many  advantages  of  multirate  signal 
processing  have  been  found,  and  techniques  have  been  applied  to  many  areas,  such  as 
telecommunications,  digital  audio  encoding/decoding,  speech  and  image  processing, 
and  geophysical  signal  processing  [Ref.  1,  2,  3]. 

Most  of  the  research  on  multirate  signal  processing,  beyond  development  and 
characterization  of  the  basic  building  blocks,  has  centered  on  filter  bank  theory  and 
multirate  Kalman  filtering  [Ref.  1,  2,  3,  4,  5,  6,  7,  8,  9,  10,  11,  12,  13,  14,  15].  A 
brief  discussion  on  each  of  these  areas  is  provided  below.  In  the  majority  of  the 
filter  bank  research,  there  is  one  input  signal  (or  vector  of  input  signals  sampled  at 
the  same  rate).  Within  the  filter  bank  structure,  the  signal  is  separated  into  different 
subbands  sampled  at  different  rates.  Signal  processing  techniques  are  applied,  and  the 
separate  signals  are  then  synthesized  into  one  output  signal  (or  vector  of  output  signals 
sampled  at  the  same  rate).  This  area  of  research  has  greatly  improved  the  knowledge 
and  understanding  of  single-input  single-output  (SISO)  and  multiple-input  multiple- 
output  (MIMO)  filter  bank  structures;  however,  it  has  not  developed  methods  for 
processing  multiple  input  signals  at  different  sampling  rates. 

In  the  area  of  statistical  multirate  signal  processing,  there  are  a  number  of 
papers  published.  One  area  that  these  papers  have  focused  on  is  the  characterization 
of  periodic  random  processes.  These  papers  present  methods  to  describe  periodic 
random  processes  in  both  the  time  and  frequency  domains  and  discuss  such  key 
concepts  as  cyclostationarity  of  periodic  random  processes.  Another  significant  area 
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of  research  in  statistical  multirate  signal  processing  is  the  development  of  state-space 
multirate  Kalman  filters.  These  Kalman  filters  have  been  applied  to  a  variety  of 
systems,  such  as  small  aircraft  orbit  control,  helicopter  passive  ranging  and  machinery 
control.  In  addition,  multirate  Kalman  filter  research  has  been  applied  to  tracking 
and  estimation  of  mobile  robots  as  well  as  recovery  of  lost  speech  packets  in  speech 
signal  processing. 

This  dissertation  approaches  statistical  multirate  signal  processing  from  the 
view  point  of  optimal  filtering  (i.  e.,  Wiener  filtering)  of  multiple  input  signals  at 
different  sampling  rates.  Following  the  investigation  of  statistical  characteristics  for 
multirate  signals,  an  optimal  estimator  for  a  desired  output  sequence  is  developed. 
This  development  is  then  generalized  to  the  linear  prediction  problem,  which  leads 
to  the  general  form  of  the  optimal  multirate  Wiener  filter.  Finally,  this  optimal  filter 
is  applied  to  the  sequential  target  classification  problem.  Much  of  the  foundation 
for  the  multirate  signal  processing  in  this  dissertation  is  derived  from  [Ref.  16,  17]. 
In  addition  the  optimal  estimation  and  prediction  problems  in  this  dissertation  are 
multirate  extensions  of  analogous  single-channel  and  multichannel  concepts  that  can 
be  found  in  [Ref.  18,  19].  The  multirate  sequential  classification  algorithm  is  derived 
from  previous  work  on  the  single-channel  and  multichannel  classification  theory  and 
algorithms  of  [Ref.  20,  21,  22], 

The  following  paragraphs  review  some  of  the  research  that  has  been  conducted 
in  multirate  signal  processing.  This  includes  multirate  theory,  multirate  filter  bank 
theory,  statistical  signal  characterization  and  multirate  Kalman  filtering.  It  is  by  no 
means  exhaustive,  and  indeed  much  work  is  still  being  conducted  today. 

A.  LITERATURE  REVIEW 

The  interest  in  multirate  signal  processing  gained  much  popularity  following 
the  1975  IEEE  Arden  House  Workshop  for  Acoustics,  Speech,  and  Signal  Processing 
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[Ref.  16].  Much  of  the  initial  research  in  multirate  signal  processing  has  focused 
on  how  to  describe  the  multirate  system,  including  the  descriptions  of  some  funda¬ 
mental  building  blocks  of  the  multirate  system  and  on  how  to  apply  multirate  signal 
processing  to  filter  bank  theory.  In  1983,  Crochiere  and  Rabiner  published  the  first 
comprehensive  book  on  multirate  signal  processing  [Ref.  16].  Their  work  presents 
the  basics  of  sampling  rate  conversion  through  decimation  and  expansion,  including 
time  and  frequency  characterization  of  multirate  signals.  The  application  that  they 
primarily  focus  upon  is  improving  efficiency  in  a  multirate  system  with  large  sampling 
rate  changes  through  multistage  structures.  These  multistage  structures  break  the 
sampling  rate  changes  into  multiple  stages,  which  can  allow  relaxation  of  filter  design 
criteria.  In  addition  to  multistage  applications,  Crochiere  and  Rabiner  also  describe 
several  forms  of  uniform  filter  banks  (e.g.,  the  discrete  Fourier  transform  (DFT)  filter 
bank  and  the  uniform  single-sideband  (SSB)  filter  bank)  where  the  signal  is  sepa¬ 
rated  into  multiple  parts  and  all  parts  are  decimated  by  the  same  value.  They  extend 
the  DFT  filter  bank  to  a  generalized  non-uniform  DFT  filter  bank,  which  allows  the 
separated  signals  to  be  decimated  by  different  values. 

Ten  years  after  the  book  by  Crochiere  and  Rabiner,  the  next  major  publication 
in  multirate  signal  processing  was  by  Vaidyanathan  in  1993  [Ref.  17].  Vaidyanathan’s 
work  addresses  the  basics  of  multirate  signal  processing  but  provides  greatly  expanded 
applications  in  filter  bank  theory.  Vaidyanathan  presents  and  thoroughly  describes 
many  types  of  multirate  filter  banks,  including  maximally  decimated  filter  banks, 
paraunitary  perfect  reconstruction  filter  banks,  linear  phase  perfectly  reconstruction 
quadrature  mirror  filter  (QMF)  banks  and  cosine  modulated  filter  banks.  His  work 
also  presents  the  concept  of  periodically  time-varying  filters,  which  plays  an  impor¬ 
tant  part  in  the  processing  of  multirate  signals  and  the  problem  of  quantization  error. 
Lastly,  Vaidyanathan  introduces  the  relationship  between  wavelet  transform  theory 
and  multirate  signal  processing  since  the  wavelet  transform  inherently  uses  nonuni¬ 
form  decimation  in  developing  its  snbbands.  Both  of  these  books  [Ref.  16,  17]  have 
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provided  the  foundation  on  which  much  of  the  research  on  multirate  signal  processing 
since  then  has  been  built. 

In  addition  to  these  two  seminal  volumes  on  multirate  signal  processing,  there 
have  been  many  contributions  in  journal  papers  and  workshop  and  conference  pub¬ 
lications.  Shenoy,  Burnside  and  Parks  extend  previous  work  on  multirate  filter  bank 
theory  to  optimum  minimax  filters  [Ref.  1] .  By  using  a  generalized  Fourier  analysis, 
they  derived  a  new  error  criteria  for  multirate  filter  design,  which  can  be  used  to 
design  the  optimum  minimax  multirate  filters  for  a  specific  input  signal  class.  In 
addition,  Chen  and  Vaidyanathan  [Ref.  4]  extend  the  concept  of  polyphase  filters  to 
describe  rational  sampling  rate  alterations.  These  polyphase  filters  are  useful  in  repre¬ 
senting  perfect  reconstruction  properties  of  multidimensional  delay-chain  systems  and 
periodicity  properties  of  decimated  periodic  signals.  Rules  for  multirate  structures 
were  also  extended  by  Evans,  Bamberger  and  McClellan  [Ref.  5]  to  multidimensional 
structures  by  Ending  greatest  common  subblattices  and  computing  coset  vectors. 

Since  the  multirate  system  contains  signals  at  different  sampling  rates,  there 
is  a  natural  periodicity  that  can  be  seen  in  the  multirate  system.  Some  basic  forms 
of  linear  periodically  time-varying  (LPTV)  filters  were  described  by  Vaidyanathan; 
however,  Saadat  Mehr  and  Chen  [Ref.  23]  present  two  methods  of  describing  LPTV 
structures,  each  of  which  consists  of  a  periodic  switch  connected  to  several  linear 
time-invariant  (LTI)  systems.  Saadat  Mehr  and  Chen  show  that  these  structures  can 
be  used  to  solve  a  general  approximation  problem  where  an  LPTV  system  with  period 
p  is  approximated  by  an  LPTV  system  with  period  p.  They  extend  their  work  on 
LPTV  structures  to  address  polyphase  and  alias-component  representations  of  LPTV 
systems  [Ref.  24] .  They  show  that  in  general  a  filter  bank  can  be  represented  by  two 
LPTV  systems  connected  in  a  cascading  manner.  They  also  show  that  a  p-channel 
filter  bank  with  period  m  can  be  represented  by  an  mp-channel  LTI  filter  bank  if  p 
and  m  are  relatively  coprime  integers. 
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Other  research  on  time-varying  systems  has  included  work  by  Phoong  and 
Vaidyanathan  [Ref.  25]  on  using  a  polyphase  approach  to  filter  banks  to  address 
unusual  properties  not  exhibited  by  LT1  filter  banks.  This  work  is  extended  in  [Ref. 
26]  to  discuss  the  factorability  of  linear  time-varying  lossless  filter  banks.  In  particu¬ 
lar,  they  show  that  all  degree-one  lossless  linear  time-varying  (LTV)  systems  can  be 
decomposed  into  a  time-dependent  unitary  matrix  and  a  lossless  dyadic-based  LTV 
system.  Vaidyanathan  expanded  his  work  on  periodic  systems  to  address  periodic  sys¬ 
tems  with  allpass  and  paraunitary  properties  [Ref.  27].  By  providing  a  state-space 
representation  of  periodic  LTI  systems  and  introducing  the  concept  of  reachability 
and  observability  in  terms  of  multirate  periodic  systems,  Vaidyanathan  shows  that 
reachability  and  observability  are  not  related  to  the  system  minimality  in  a  simple 
way,  unlike  traditional  state-space  linear  systems.  In  [Ref.  28],  Gadre  and  Patney 
address  issues  associated  with  aliasing  cancelation  and  perfect  reconstruction  within 
a  vector  context,  and  they  define  a  vector  multirate  system.  These  definitions  lead 
to  conditions  whereby  a  vector  LPTV  system  can  become  a  time-invariant  system. 
Further,  Spurbeck  and  Scharf  [Ref.  29]  applied  spectral  factorization  techniques  to 
the  filter  design  for  periodically  correlated  time  series. 

Another  approach  to  representing  periodic  systems  was  introduced  by  Misra 
[Ref.  30].  Misra  shows  that  many  results  from  linear  time  invariant  theory  can  be 
extended  to  periodic  systems.  A  necessary  condition  for  these  results  is  that  an 
equivalent  time  invariant  system  must  be  found  for  the  periodic  system.  Misra  pro¬ 
vides  a  numerically  reliable  and  simple  procedure  to  find  the  equivalent  time  invariant 
systems,  and  he  shows  that  a  minimal-order  generalized  state-space  description  can 
always  be  found. 

The  research  cited  up  to  this  point  has  focused  on  deterministic  signal  pro¬ 
cessing.  In  addition,  the  primary  focus  has  centered  on  filter  bank  theory  and  how 
to  best  separate  a  signal  into  multiple  parts  that  are  sampled  at  different  rates  per¬ 
form  any  necessary  signal  processing  and  then  resynthesizing  the  processed  signals. 
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In  many  signal  processing  techniques,  however,  applications  involve  signals  that  are 
random  processes.  Consequently,  some  work  has  been  done  to  appropriately  describe 
the  statistical  characterizations  of  multirate  signal  processing  techniques. 

Since  multirate  systems  are  inherently  periodic  in  nature,  signals  within  the 
multirate  system  are  usually  not  wide-sense  stationary  (WSS).  However,  under  ap¬ 
propriate  conditions,  stationary  input  signals  can  exhibit  cyclostationary  or  wide- 
sense  cyclostationary  (WSCS)  characteristics.  In  [Ref.  31]  and  [Ref.  32],  Gardner 
presents  methods  on  how  to  spectrally  characterize  Nth  order  cyclostationary  signals 
and  how  to  exploit  spectral  redundancies  that  might  exist  in  cyclostationary  signals. 
Related  work  in  characterization  of  LTV  systems  was  conducted  by  Akkarakaran  and 
Vaidyanathan  [Ref.  33]  through  the  use  of  bifrequency  and  bispectrum  maps.  The 
bifrequency  map  is  a  two-dimensional  Fourier  transform  used  to  characterize  the  LTV 
system  while  the  bispectrum  map  is  a  two-dimensional  Fourier  transform  that  charac¬ 
terizes  non-stationary  random  processes.  In  [Ref.  34]  and  [Ref.  35] ,  Therrien  presents 
methods  for  defining  correlations  functions  and  power  spectra  for  multirate  signals. 
This  is  extended  by  Therrien  to  non-stationary  random  processes  in  [Ref.  36]  and  in 
[Ref.  37], 

Research  into  cyclic  higher-order  statistics  of  multirate  signals  has  been  con¬ 
ducted  by  Napolitano  [Ref.  38].  Napolitano  uses  cyclic  higher-order  statistics  to 
derive  the  input-output  relations  for  MIMO  linear  almost-periodically  time-variant 
systems  excited  by  cyclostationary  inputs.  In  addition,  he  presents  a  sufficient  condi¬ 
tion  on  the  sampling  rate  to  prevent  aliasing  when  reconstructing  cyclic  higher-order 
statistics  of  continuous  signals  from  sampled  signals.  Research  into  the  problem  of 
avoiding  aliasing  in  the  cyclic  higher  order  spectra  of  a  decimated  time  series  is  ex¬ 
tended  by  Izzo  and  Napolitano  in  [Ref.  39]  to  a  generalized  form  that  avoids  aliasing 
and  imaging  effects  of  a  time  series  decimated  by  a  fractional  factor. 

LI sing  cyclostationary  spectral  analysis,  Ohno  and  Sakai  [Ref.  6]  developed  a 
method  of  deriving  multirate  optimal  biorthogonal  FIR  filter  banks  that  minimizes 
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the  time-averaged  mean-square  error  (TAMSE)  after  the  high-frequency  subband  is 
removed.  In  order  to  properly  describe  the  filter  bank,  cyclostationary  spectral  anal¬ 
ysis  is  used  since  the  output  of  the  filter  bank  is  cyclostationary  for  a  wide-sense 
stationary  input. 

Other  statistical  methods  are  used  by  Yurdakul  and  Dundar  [Ref.  40]  to 
estimate  the  quantization  error  of  FIR-based  multirate  systems.  This  is  practical  for 
implementation  of  multirate  filters  since  the  quantization  of  filter  coefficients  leads  to 
quantization  errors  in  the  output  signal  of  any  discrete-time  system. 

One  technique  of  statistical  multirate  signal  processing  researched  by  Sathe 
and  Vaidayanathan  [Ref.  41]  is  adaptive  filtering.  They  use  the  statistical  proper¬ 
ties  of  signals  in  LTV  systems  to  develop  an  adaptive  filter  structure  that  is  useful 
for  identification  of  band-limited  channels.  A  matrix  form  of  this  adaptive  filter  is 
shown  to  provide  better  performance  in  terms  of  lower  error  energy  than  a  traditional 
adaptive  filter  but  suffers  from  high  computational  burden. 

Outside  of  the  signal  processing  literature,  the  most  prominent  technique  in 
statistical  multirate  signal  processing  available  is  Kalman  filtering.  Many  applications 
of  and  techniques  for  implementing  multirate  Kalman  Liters  have  been  investigated. 
Early  work  on  multirate  Kalman  filtering  by  Andrisani  and  Gau  [Ref.  7]  uses  two 
different  Kalman  Liters  in  parallel.  One  of  the  Kalman  Liters  processes  the  high  rate 
measurement,  at  a  reduced  order.  The  second  Kalman  Liter  processes  the  residuals 
of  the  Lrst  Kalman  Liter  in  conjunction  with  the  low-rate  measurement.  Andrisani 
and  Gau  were  able  to  reduce  the  computational  complexity  of  their  multirate  Kalman 
Liter  by  developing  a  suboptimal  version,  although  a  slight  performance  penalty  was 
incurred. 

Much  work  in  Kalman  Lltering  has  been  advanced  by  Bor-Sen  Chen,  along 
with  You-Lin  Chen  and  Chin- Wei  Lin  [Ref.  2,  3,  8,  9].  Through  their  work,  the 
multirate  Kalman  Liter  has  been  applied  to  modeling  of  autoregressive  (AR)  and 
moving-average  (MA)  models  through  interpolation  and  estimation  of  the  values  of 
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the  AR  or  MA  stochastic  signals.  In  addition,  these  researchers  have  developed  meth¬ 
ods  for  optimal  signal  reconstruction  in  noisy  filter  banks  and  developed  a  multirate 
Kalman  filter  that  models  system  and  channel  noise.  In  addition  to  using  a  multirate 
Kalman  reconstruction  filter,  they  derived  a  sample-interpolation  algorithm  useful  for 
the  recovery  of  missing  speech  packets. 

Further  work  on  Kalman  filtering  was  conducted  by  De  Leon,  Kober,  Krumvieda 
and  Thomas  [Ref.  10].  By  separating  signals  into  subbands  and  then  implementing 
multiple  Kalman  Liters,  a  target  can  be  tracked  while  reducing  the  computational  rate 
and  sensitivity  to  noisy  measurements.  Ni  et  al.  [Ref.  11]  also  developed  multirate 
Kalman  Liters  that  are  less  sensitive  to  noise.  Their  model  combines  the  statistical 
model  of  the  input  signal  with  a  multichannel  representation  of  the  subband  signal. 
This  form  of  the  multirate  Kalman  Liter  provides  the  minimum  variance  reconstruc¬ 
tion  of  the  input  signal,  provided  that  the  input  signal  is  embedded  in  the  state  vector. 
Other  applications  of  the  multirate  Kalman  Liter  include  orbit  control  of  small  air¬ 
craft  [Ref.  13],  helicopter  passive  ranging  [Ref.  42,  43]  and  C.N.C.  machining  control 
[Ref.  12].  In  addition,  Tornero  et  al.  [Ref.  14]  combine  the  multirate  Kalman  Li¬ 
ter  with  a  multirate  LQG  controller  to  be  used  in  self-location  and  path-tracking  in 
mobile  robots  applications. 

Another  signiLcant  area  of  research  employing  multirate  Kalman  Lltering  is 
the  multiresolution  multirate  (MRMR)  estimation  problem.  One  of  the  fundamental 
tools  of  MRMR  techniques  is  the  wavelet  transform.  The  goal  of  MRMR  estimation 
is  to  improve  the  overall  estimation  capability  by  exploiting  features  at  different 
resolutions  and  sampling  rates.  Early  work  in  this  area  was  conducted  by  Basseville 
et  al.  [Ref.  44]  on  modeling  and  estimation  of  multiresolution  statistical  processes. 
Chou  et  al.  [Ref.  45,  46]  extended  this  work  to  recursive  estimation  and  the  kalman 
Liter.  More  recently,  Cristi  and  Tummala  [Ref.  15]  extended  the  work  even  further 
by  developing  a  recursive  MRMR  Kalman  Liter. 

In  addition  to  multirate  Kalman  Lltering,  some  research  on  multirate  estima- 


tion  is  available.  Haddad,  Bernstein  and  Huang  developed  reduced  order  estimators 
[Ref.  47].  These  estimators  were  derived  using  a  system  of  equations  consisting  of  one 
modified  Riccati  equations  and  two  modified  Lyapunov  equations,  and  each  system 
of  equation  corresponds  to  a  subinterval  of  the  period  associated  with  the  multirate 
system.  Another  method  of  estimation  was  introduced  by  Hong  [Ref.  48,  49].  His 
method  focuses  on  using  high-rate  measurements  to  make  low-rate  estimations.  By 
using  filter  banks  and  taking  advantage  of  the  lowpass  filtering  effect  of  the  filter  bank, 
improved  approximations  of  the  low  rate  estimation  were  obtained.  In  addition,  a 
good  approximation  of  the  original  measurements  was  obtained  through  orthogonal 
transformation  by  using  highpass  filtering  and  lowpass  filtering.  This  thesis  had  its 
beginning  in  early  work  by  Therrien  and  his  students  and  there  have  been  several 
publications  along  the  way.  In  addition,  some  of  the  early  work  discussed  here  pro¬ 
vided  the  groundwork  for  a  companion  thesis  [Ref.  50]  that  extended  the  work  along 
various  other  directions,  most  particularly  to  two-dimensional  signals  with  applica¬ 
tions  to  super-resolution  image  reconstruction.  In  [Ref.  51],  Cristi,  Koupatsiaris  and 
Therrien  present  a  method  of  multirate  estimation  for  a  two-signal  input  along  with 
a  quantitative  analysis  of  reduction  in  mean-square  error.  The  multirate  estimator 
is  extended  to  m-signals  by  Kuchler  and  Therrien  in  [Ref.  52]  (preliminary  work  for 
this  dissertation),  and  a  recursive  means  of  finding  the  filter  coefficients  and  error 
variance  as  signals  are  added  to  the  multirate  system  is  presented.  In  addition,  Ther¬ 
rien  presents  a  method  of  linear  prediction  for  multirate  systems  in  [Ref.  53]  along 
with  a  Levinson-type  recursion  for  finding  the  prediction  parameters  (in  conjunction 
with  this  dissertation).  Further  work  by  Therrien  and  Hawes  in  [Ref.  54]  and  [Ref. 
55]  present  a  method  of  least  mean  squares  (LMS)  calculation  for  multirate  systems 
with  two  input  signals.  Lastly,  some  techniques  used  in  the  one-dimensional  multirate 
system  have  been  extended  to  the  two-dimensional  multirate  system  by  Scrofani 
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and  Therrien  in  [Ref.  56]  and  [Ref.  57].  They  use  multirate  estimation  algorithms 
to  reconstruct  a  high-resolution  image  from  a  set  of  low  resolution  images  that  are 
subpixel  translated  (super-resolution). 

B.  DISSERTATION  OUTLINE 

From  the  literature  review,  it  is  apparent  that  the  development  of  multirate 
filter  banks  has  matured  greatly.  These  filter  banks,  however,  apply  to  SISO  systems 
or  MIMO  systems  where  the  sampling  rates  of  all  the  input  and  output  signals  are 
the  same.  This  dissertation  begins  by  developing  a  general  framework  for  multirate 
signal  processing.  The  characterization  of  random  processes  as  they  are  applied  to 
decimators,  expanders  and  filters  are  presented.  Some  LPTV  filter  bank  realizations 
are  discussed  and,  in  particular,  an  LPTV  filter  bank  that  allows  the  input  and  output 
to  be  observed  at  different  sampling  rates  is  presented.  Further  simplification  of  signal 
presentation  is  presented  through  the  use  of  matrix  forms  and  Kronecker  products. 
Whereas  the  majority  of  statistical  multirate  signal  processing  has  focused  on  the 
Kalman  filter,  this  dissertation  approaches  statistical  multirate  signal  processing  from 
the  point  of  view  of  the  more  general  multirate  Wiener  filter  and  its  extensions. 

The  remainder  of  the  dissertation  is  organized  as  follows.  In  Chapter  II,  the 
basic  framework  for  describing  a  multirate  system  is  developed.  This  includes  de¬ 
scribing  the  appropriate  indexing  schemes  necessary  to  adequately  describe  different 
signals  that  are  observed  at  different  sampling  rates  and  how  they  can  be  referenced 
to  each  other.  This  chapter  also  describes  the  different  components  available  to 
multirate  signal  processing  systems.  In  particular,  the  expander,  decimator,  linear 
time-invariant  (LTI)  filter  and  linear  periodically  time-varying  (LPTV)  filter  are  de¬ 
scribed.  Further,  the  correlation  equations  for  these  components  are  presented.  In 
Chapter  III,  the  multirate  Wiener-Hopf  equations  for  optimal  filtering  are  developed. 
The  multirate  optimal  filter  is  presented  both  in  its  direct  form  and  in  its  innovations 
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representation  form.  Simulations  are  conducted  to  analyze  the  performance  of  the 
multirate  optimal  filter.  In  Chapter  IV,  the  multirate  optimal  filter  is  focused  on  the 
problem  of  linear  prediction.  A  convenient  formulation  of  this  problem  is  presented 
and  a  novel  efficient  algorithm  for  solving  the  multirate  linear  prediction  coefficients  is 
developed.  In  Chapter  V,  the  signal  processing  application  of  sequential  classification 
is  developed  for  a  multirate  system.  Simulations  are  conducted  to  explore  factors 
that  affect  the  overall  capabilities  of  multirate  sequential  classification.  In  addition, 
performance  results  for  a  single  channel,  a  multichannel  and  a  multirate  system  are 
compared.  Finally,  in  Chapter  VI,  the  research  and  conclusions  of  this  dissertation 
are  summarized  and  areas  of  further  research  are  suggested. 
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II.  ANALYSIS  OF  MULTIRATE  SIGNALS  AND 

SYSTEMS 


A.  GENERAL  CONCEPTS 

Signal  processing  techniques  have  been  well  developed  for  processing  a  single- 
input  and  single-output  system,  or  for  processing  a  multiple-input  and  multiple- 
output  system  where  all  the  signals  are  received  at  the  same  sampling  rate;  the  latter 
is  referred  to  as  multichannel  signal  processing.  In  some  signal  processing  situations, 
however,  observed  signals  are  sampled  at  various  different  rates,  and  it  may  be  desired 
to  process  or  exploit  these  signals  together  to  perform  various  operations,  such  as  esti¬ 
mation,  prediction,  detection  or  classification.  To  provide  a  basis  for  addressing  some 
of  these  problems,  some  new  theory  needs  to  be  developed  that  extends  the  theory  for 
single-rate  signals  and  systems.  The  theory  developed  for  multirate  statistical  signal 
processing  should  reduce  to  the  single-channel  and  single-rate  multichannel  problems 
as  special  cases.  In  this  chapter,  some  basic  concepts,  notation,  and  terminology 
for  describing  stochastic  multichannel  systems  and  signals  are  introduced.  In  other 
words,  this  chapter  presents  the  mathematical  tools  to  be  used  in  the  remainder  of 
the  research. 

A  representation  of  a  general  multirate  system  is  shown  in  Fig.  2.1.  From  Fig. 
2.1  it  can  be  seen  that  each  signal  has  an  associated  sampling  rate,  Ft  (Hz),  where 
the  letter  i  represents  a  particular  signal.  These  sampling  rates  are  generally  differ¬ 
ent,  although  some  may  be  identical.  Since  these  signals,  in  general,  have  different 
sampling  rates,  their  observations  do  not  necessarily  occur  at  the  same  point  in  con¬ 
tinuous  time  and  their  indices  m,  may  not  be  aligned.  For  example,  suppose  a  signal 
x  is  sampled  at  3000  Hz  and  a  signal  y  is  sampled  at  2000  Hz.  If  both  signals  have 
a  sample  that  occurs  at  t  —  0,  then  the  next  observation  of  the  signal  x  would  occur 
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Figure  2.1.  Multirate  System  Concept. 


at  0.33  milliseconds  and  the  next  observation  of  y  would  occur  at  0.5  milliseconds. 
Observations  of  both  x  and  y  would  not  coincide  until  1  millisecond  after  the  first 
observation  of  both  signals,  this  concept  is  illustrated  in  Fig.  2.2.  A  diagram  such  as 
that  provided  in  Fig.  2.2  will  be  called  a  sampling  pattern. 


Fx  =  3000  Hz 

x  •  •  •  • . . . 

— f - 1 - 1 - 1 - 1 - 1 - r*-  time 

0  0.5  ms  1  ms 

y  •  •  «... 

F  =  2000  Hz 

Figure  2.2.  Effect  of  Different  Sampling  Rates  on  Observation  Times. 

In  order  to  provide  a  common  framework  to  which  all  the  signals  can  be  refer¬ 
enced,  to  account  for  the  effects  of  sampling  signals  at  different  rates,  a  structure  is 
defined  here  that  will  be  called  the  fundamental  layer.  This  structure  does  not  neces¬ 
sarily  correspond  to  any  physical  portion  of  the  multirate  system.  It  is  a  conceptual 
tool  that  is  used  to  allow  common  referencing  of  all  the  signals  within  the  multirate 
system. 
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A  key  element  of  the  fundamental  layer  is  that  it  has  a  sampling  rate  that  is 
a  multiple  of  the  sampling  rates  for  all  signals  in  the  system.  This  allows  all  of  the 
signals  to  be  described  explicitly  in  terms  of  fundamental  layer  parameters.  Return¬ 
ing  to  the  example  of  two  observed  signals  x  and  y  with  sampling  rates  Fx  =  3000 
Hz  and  Fy  =  2000  Hz,  respectively,  the  concept  of  the  fundamental  sampling  rate 
is  illustrated.  In  order  to  represent  both  signals  in  a  common  framework  (the  fun¬ 
damental  layer)  for  analysis,  each  signal  is  represented  at  a  higher  rate  which  is  an 
integer  multiple  of  the  actual  sampling  rate.  For  efficiency,  the  lowest  such  rate  that 
can  be  applied  to  all  signals  is  used.  In  the  example  just  described,  this  sampling  rate 
would  be  6000  Hz,  and  is  illustrated  by  the  tick  marks  on  the  time  axis  of  Fig.  2.2. 
This  sampling  rate  will  be  called  the  fundamental  rate ,  F  and  its  inverse  T  =  l/F 
will  be  called  the  system  clock  rate.  The  fundamental  rate  is  the  minimum  sampling 
rate  necessary  to  describe  all  signals  in  the  multirate  system.  It  is  assumed  that  each 
sampling  rate  Ft  is  integer- valued;  thus  the  following  definition  can  be  made. 


Definition  1.  The  Fundamental  Rate,  F ,  is  given  by 

~F  —  LCM(Fi,  F2, . . .),  (2.1) 

where  LCM(  )  denotes  the  least  common  multiple  and  F1;  F2,  . . .  represent  the  sam¬ 
pling  rates  of  all  the  signals  in  the  system. 

Each  actual  signal  in  the  system  can  be  associated  with  (a  possibly  fictitious) 
signal  in  the  fundamental  layer  which  is  sampled  at  the  fundamental  rate  F.  Thus 
each  signal  in  the  system  can  be  thought  of  as  a  down-sampled  or  decimated  version 
of  some  signal  at  the  fundamental  rate.  The  decimation  factor  associated  with  the 
ith  signal  is  given  by 

Ki  =  T  (2.2) 

F  i 

i.e.,  it  is  the  ratio  of  the  fundamental  sampling  rate  to  the  signal  sampling  rate  for 
the  ith  signal  of  interest.  While  notation  such  as  x \mx  is  used  to  describe  a  signal  at 
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some  sampling  rate  Fx,  the  notation  x[n]  will  be  used  to  describe  the  corresponding 
signal  in  the  fundamental  layer.  Thus  these  two  signals  are  related  by 


x[mx 


x[n 


n=Kxmx 


x[Kxmx 


For  a  system  conceptualized  as  in  Fig.  2.1,  a  global  periodicity  exists.  The 
system  period  represents  the  minimum  number  of  time  steps  at  the  fundamental  rate 
necessary  to  establish  a  repetitive  sampling  pattern  for  all  signals  in  the  multirate 
system. 


Definition  2.  The  System  Period,  K,  is  the  minimum  common  digital  period  of  all 
elements  of  the  digital  system.  The  system  period  is  given  by 


K  =  LCM(Ad,  K2, . . .), 


(2.3) 


where  K\,  K2,  ■  ■  ■  represent  the  decimation  factors  of  all  the  signals  in  the  system. 

For  the  ith  signal,  the  number  of  samples  per  period  is  given  by 

K  , 

Mi  =  w.  (2,4) 

-K-i 

While  K  represents  the  number  of  time  steps  at  the  system  rate  corresponding  to  one 
period,  M*  represents  the  number  of  time  steps  at  the  (ith)  signal  rate  corresponding 
to  one  period. 

In  order  to  demonstrate  the  concepts  introduced  thus  far,  the  following  exam¬ 
ple  of  a  two-channel  multirate  system  is  provided.  Consider  a  multirate  system  that 
has  two  signals  x  and  y  with  sampling  rates  Fx  =  3000  Hz  and  Fy  =  2000  Hz,  respec¬ 
tively.  From  (2.1),  the  fundamental  rate  for  this  system  is  F  —  LCM(3000,  2000)  = 
6000  Hz.  Then,  from  (2.2),  the  signal  decimation  factors  are  calculated  to  be  Kx  =  2 
and  Ky  =  3.  Finally,  from  (2.3),  the  system  period  is  found  to  be  K  =  6.  Figure 
2.3  shows  the  relationships  among  these  terms.  It  can  be  seen  from  Fig.  2.3  that  the 
decimation  factors  Kx  and  Ky  represent  the  number  of  time  steps,  in  the  fundamental 
layer,  between  successive  observations  (samples)  of  a  particular  signal.  Notice  that 
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K  =  6  K 


y  O  .  •  Q  .  .  Q — — • — O— - O — «■ — «— O — • 


Figure  2.3.  Illustration  of  Decimation  Factors  and  System  Period. 

the  sampling  pattern  repeats  every  K  =  6  time  steps  at  the  fundamental  rate.  This 
also  corresponds  to  Mx  =  3  samples  at  the  rate  of  x  and  My  =  2  samples  at  the  rate 
of  y.  Thus  given  a  time  n,  the  relative  position  within  the  sampling  pattern  is  the 
same  for  all  times  n  +  iK  where  i  is  an  integer.  This  relative  position  will  be  called 
the  system  phase  of  the  system  (not  to  be  confused  with  the  phase  of  its  Fourier 
transform).  The  system  phase,  denoted  by  k  can  be  determined  for  any  time  n  by 
writing 

n  =  mK  +  k, 

where  m  =  \n/K\.  Thus  the  following  definition  can  be  made. 

Definition  3.  The  System  Phase,  k,  is  given  by 

n  =  k  (mod  K),  (2.5) 

where  n  is  associated  time  at  the  fundamental  layer,  and  k  is  the  residual  modulo  K . 


B.  MULTIRATE  SIGNAL  PROCESSING 

Some  common  forms  of  discrete-time  signals  can  be  directly  related  to  an 
associated  continuous-time  signal.  This  will  be  shown  using  the  sampling  rates  of 
the  observed  signals  of  a  multirate  system  and  the  definition  of  the  fundamental  rate. 
Given  a  continuous  time  signal  x(t),  the  discrete-time  signal  for  a  signal  sampled  at 
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Figure  2.4.  Continuous  and  Multirate  Signals:  (a)  continuous  signal,  (b)  discrete 
signal  at  fundamental  sampling  rate,  (c)  discrete  signal  at  observed  sampling  rate. 

the  rate  Fx.  would  be  Xi[mXi ]  =  x(mXiTXi ),  where  TXi  =  — —  and  i  —  1, . . . ,  N,  where 

I’xi  _ 

N  is  the  number  of  signals  in  the  system.  Calculating  the  fundamental  rate  as  F  = 
LCM(FXl,  FX2 , . . . ,  FXn ) ,  the  discrete-time  signal  associated  with  the  fundamental  rate 
would  be  x[n]  =  x(nT),  where  T  —  =  is  the  system  clock  rate.  All  signals  are 
referenced  to  the  initial  time  of  t  —  0,  thus  a;[0]  =  x(0Tx)  =  x(0).  For  clarity, 
parentheses  are  used  to  indicate  the  time  dependence  for  a  continuous  signal,  and 
brackets  are  used  to  indicate  the  time  step  for  a  discrete-time  signal. 

The  relationships  among  some  various  signals  are  depicted  in  Fig.  2.4.  The 
original  continuous  signal  is  denoted  by  x{t).  The  signal  in  the  fundamental  layer  is 
x [n] ,  which  may  not  exist  in  reality,  and  the  signal  which  is  sampled  at  the  given  rate 
FXl  is  then  x\ [mXl],  which  is  a  decimated  version  of  x[n\. 
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Some  common  multirate  signal  processing  operations  are  now  discussed  in 


detail. 


1.  Decimation 


The  concept  of  decimation,  or  downsampling,  was  introduced  earlier  in  Section 
II.A  in  conjunction  with  the  definition  of  signals  in  the  fundamental  layer.  In  this 
section,  decimation  is  discussed  in  a  more  general  sense,  as  it  may  relate  to  various 
signals  that  may  exist  in  a  multirate  system.  Here  decimation  is  discussed  from  the 
point  of  view  of  obtaining  one  signal  within  the  system  from  another. 

The  process  of  decimation  eliminates  data  points  in  a  signal  vector  in  order  to 
reduce  the  sampling  rate.  The  number  of  data  points  removed  is  dependent  on  the 
decimation  factor.  Given  an  integer-valued  decimation  factor,  L,  only  one  of  every  L 
consecutive  samples  is  retained. 


y[my\  =  x[Lmy]  ,  for  rny  =  . . . ,  —1,  0, 1, ... . 


(2.6) 


The  implied  sampling  rate  for  y  is  thus  also  reduced  by  a  factor  of  L,  i.e. ,  Fy  =  Fx/L. 
The  process  of  decimation  is  represented  in  block  diagram  form  in  Fig.  2.5,  where  L 
represents  the  decimation  factor  and  my  is  the  index  of  the  decimated  signal,  such 
that  rriy  is  the  set  of  integers  Z  =  {. . . ,  —1,  0, 1, . . .}.  The  associated  index  of  the 
input  signal  is  mx  =  Lmy. 


*  [  mx  ] - »  i  l  - *y  \my  ] 


Figure  2.5.  Decimation. 


It  is  well  known,  [Ref.  17,  58],  that  the  frequency  spectrum  Y( coy)  of  the 
decimated  signal  is  related  to  the  frequency  spectrum  X(ux)  of  the  original  signal  by 
the  formula 


(2.7) 


The  effect  of  decimation,  in  the  frequency  domain,  is  to  stretch  the  frequency  band 
as  shown  in  Fig.  2.6.  Notice  that  if  the  spectrum  of  the  original  signal  is  nonzero  for 
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|u;|  >  — ,  aliasing  will  occur  due  to  undersampling.  In  most  cases  this  is  undesirable, 
L 

although  in  some  cases  aliasing  can  be  exploited  to  obtain  useful  results.  To  eliminate 
aliasing,  the  data  sequence  may  be  pre-hltered  with  an  (ideal)  low  pass  filter  with  a 

7T 

cutoff  frequency  at  |u;|  =  — . 


Figure  2.6.  Frequency  Spectrum  for  Decimation 

Notice  that  given  a  signal  x\n]  at  the  fundamental  rate,  it  is  possible  to  generate 
L  different  decimated  signals  by  applying  a  time  shift  to  the  signal  at  the  fundamental 
rate  before  decimation.  These  L  signals  all  have  the  same  rate  and  are  defined  by 


y^[m]  =  x[Lmy  +  /],  l  —  0, 1, . . . ,  L  —  1. 
The  generation  of  these  signals  is  illustrated  in  Fig.  2.7, 


rM 


lL 


JO 


M 


Figure  2.7.  Decimation  with  Time  Shift  l. 


where  zl  represents  a  time  shift  of  /  units  at  the  input  rate.  When  all  of  the  L  possible 
decimation  signals  are  formed,  this  process  is  known  as  complete  decimation.  Figure 
2.8  illustrates  complete  decimation  for  a  signal  x[n]  at  the  fundamental  rate. 

2.  Expansion 

The  opposite  of  decimation  is  the  process  of  expansion.  Expansion,  or  up- 
sampling,  increases  the  sampling  rate  by  inserting  zeros  between  the  data  points  of  a 
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Figure  2.8.  Complete  decimation  illustrated  for  L  =  3:  (a)  discrete  signal  at  funda¬ 
mental  sampling  rate,  (b)  discrete  signal  at  decimated  sampling  rate,  no  unit  shift 
(1  =  0),  discrete  signal  at  decimated  sampling  rate,  one  unit  shift  (7  =  1),  (d)  discrete 
signal  at  decimated  sampling  rate,  two  unit  shifts  (1  =  2). 
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signal.  The  sampling  rate  is  increased  by  a  factor  of  /  if  /  —  1  zeros  are  inserted  after 
every  data  point.  This  process  is  defined  by 


y[my] 


I  x[my/I\  for  my  div  / 
I  0  otherwise 


(2.8) 


where  the  expression  amy  div  J”  signifies  that  my  is  divisible  by  I,  i.e. ,  my/I  is  an 
integer.  The  process  is  represented  in  block  diagram  form  in  Fig.  2.9.  The  frequency 


*K] 


t  / 


Figure  2.9.  Expansion. 


spectrum  of  the  expanded  signal  is  related  to  the  original  signal  by 


1>.)  =  X(IUy). 


(2.9) 


The  effect  of  expansion  in  the  frequency  domain  is  to  compress  the  frequency  band 
as  shown  in  Fig.  2.10.  This  compression  brings  additional  copies  of  the  original 


X(^) 


Figure  2.10.  Frequency  spectrum  for  Expansion. 


spectrum,  known  as  “images”,  into  the  frequency  band  of  interest  ( — 7r,  7t).  In  order 
to  produce  a  rate  change  without  distortion  the  expander  must  be  followed  by  an 

7 T 

(ideal)  low-pass  filter  with  cutoff  frequency  at  \cu\  —  In  the  time  domain,  filtering 
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after  upsampling  has  the  effect  of  “filling  in”  the  zeros  with  “interpolated”  values  of 
the  signal. 


3.  Rate  Changes 


Any  rate  change  by  a  rational  factor,  I/L,  can  be  achieved  by  expansion, 


low-pass  filtering  and  then  decimation  of  a  signal.  (It  is  assumed  that  /  and  L 


are  coprime,  i.e.,  these  integers  have  no  common  factors.)  The  structure  for  this 
rate  change  system  is  shown  in  Fig.  2.11.  The  input  x[mx]  is  at  a  rate  Fx.  After 
expansion,  the  rate  changes  to  FXI .  The  low-pass  filter  which  is  designed  to  operate 
at  this  new  rate  ( FXI )  then  serves  both  to  eliminate  the  images  due  to  expansion  and 
provide  bandlimiting  so  the  decimated  signal  will  not  be  aliased.  In  order  to  avoid 
aliasing,  however,  the  cutoff  frequency  of  the  low  pass  filter  must  be  less  than  or  equal 


Fy  —  Fx-  I/L. 


Observe  that  the  system  rate  for  this  simple  system  is  F  =  FXI  and  that 
the  decimation  factors  are  given  by  Kx  =  /  and  Ky  =  L.  Thus  a  rate  change  is 
accomplished  by  expanding  up  to  the  system  rate,  filtering  and  then  decimating. 


rate  Fx 


rate  FI  rate  Fv  =  FT  /  L 

x  y  x 


Figure  2.11.  Rate  Change  by  a  Rational  Factor  I/L. 
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4.  Filtering 

Within  the  framework  of  multirate  signal  processing,  there  are  two  broad  cat¬ 
egories  of  filters  that  are  most  useful.  These  filters  are  1)  linear  time-invariant  filters 
and  2)  linear  periodically  time-varying  filters,  and  are  discussed  in  the  following  sub¬ 
subsections. 

a.  Linear  Time-Invariant  Filters 

Linear  time-invariant  Liters  process  a  signal  at  some  sampling  rate  Fx 
and  produce  an  output  at  the  same  rate.  The  input  x  and  output  y  are  related  by 
convolution  as 


OO  OO 

y[m]  —  h[k]x[m  —  k]  —  x[m]h[n  —  m }  (2-10) 

k=— oo  m=— oo 

where  h[m]  is  called  the  impulse  response  and  has  the  same  rate  as  x  and  y.  The 
system  is  causal  if  an  only  if  h[m]  —  0  for  rn  <  0.  In  this  case  the  lower  limit  of  the 
variable  m  in  (2.10)  can  be  changed  from  —  oo  to  0. 

The  system  is  stable  if  and  only  if 


OO 

m=— oo 


<  (X). 


The  system  is  said  to  have  finite  impulse  response  (FIR)  if  h[m]  is  a  finite-length 
sequence  and  said  to  have  infinite  impulse  response  (IIR)  otherwise.  In  the  case  of 
an  FIR  Liter  (2.10)  can  be  reduced  to  a  Lnite  sum  and  the  system  is  always  stable. 
b.  Linear  Periodically  Time- Varying  Filters 
Linear  periodically  time- varying  (LPTV)  Liters  are  important  in  multi¬ 
rate  applications  because  of  the  inherent  periodic  nature  of  the  multirate  system.  For 
example,  in  optimum  Lltering  applications  (see  Chapter  IV),  the  system  periodicity 
requires  that  the  Liter  coefficients  change  in  a  periodic  fashion.  This  is  illustrated  in 
Fig.  2.12  where  the  Liter  output  y[my\  is  intended  to  estimate  another  signal  d[my]. 
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Figure  2.12.  LPTV  Filter 


The  sampling  patterns  for  x,  y  and  d  are  shown  in  Fig.  2.12  together 
with  a  time  scale  labeled  ‘n’  that  represents  the  system  rate  in  the  fundamental  layer. 
The  filtering  is  assumed  to  be  causal  so  that  at  any  selected  time  step  only  values 
of  the  input  occurring  at  the  same  time  or  earlier  on  the  system  time  scale  are  used. 
From  Fig.  2.12  it  can  be  seen  that  at  time  my  =  6  the  input  sequence  x  and  the  output 
sequence  y  are  aligned,  however  at  time  my  =  7  the  input  sequence  lags  the  output 
sequence.  The  difference  in  system  clock  time  exhibited  at  these  points  implies  there 
will  be  a  difference  in  correlation  between  input  and  output.  Therefore,  to  estimate 
the  output  properly  at  time  my ,  a  new  set  of  filter  coefficients  different  from  those 
at  time  my  =  7  are  needed,  although  the  filter  coefficients  are  operating  on  the  same 
input  sequence. 

Consider  the  following  specific  example.  Given  a  causal  LPTV  filter 
with  a  Liter  order  of  four,  the  optimal  estimate  of  d  at  time  my  =  6  is 

y[  6]  =  /z(0)[OM4]  +  M°)[1M3]  +  /r(0)[  2]x[2]  +  h^[3]x[l],  (2.11) 

where  h^ll],  l  =  0, 1,  2,  3  are  the  optimal  Liter  coefficients  for  estimating  d[ 6].  At  the 
next  desired  estimate  time,  rny  =  7,  no  new  input  data  has  been  observed.  Therefore, 
the  estimate  y[ 7]  also  is 

y[  7]  =  /r(2)[  0]x[4]  +  h<-2)[l]x[3]  +  h(2)[  2]x[2]  +  h^[3]x[l],  (2.12) 
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but  h^[l],  l  =  0, 1,  2,  3  is  a  different  set  of  coefficients  since  the  filter  is  estimating 
d{ 7],  instead  of  d[6]. 

At  time  my  =  8,  a  new  observation  of  the  input  has  already  occurred. 
Therefore,  the  filter  equation  at  this  time  instance  is 

y[8]  =  [0]x[5]  +  h(4)[  l]x[4]  +  hw[2]x[3}  +  hw[3]x[2].  (2.13) 

Again  the  coefficients  are  different  from  those  occurring  in  (2.11)  and  (2.12).  Finally, 
at  time  rny  =  9,  the  input  observation  occurs  at  the  same  time  as  the  output  esti¬ 
mate  is  calculated.  This  is  the  same  system  phase  at  the  time  instance  of  my  =  6. 
Therefore,  the  optimal  filter  coefficients  at  time  my  =  9  are  the  same  as  those  at  time 
my  =  6. 

Although  this  example  uses  specific  time  steps  to  show  how  an  LPTV 
filter  functions,  a  generalized  filter  equation  can  be  found.  Let  x  and  y  be  the  repre¬ 
sentations  of  x  and  y  in  the  fundamental  layer,  and  let  [i]  be  a  set  of  time- varying 
filter  coefficients.  The  most  generic  linear  filter  can  be  written  as 

OO 

y[n\  —  ^  h^[i]x[n  —  *].  (2-14) 

i=— oo 

Observations  of  x  only  occur  when  its  argument  is  a  multiple  of  Kx  however.  Thus 
the  right-hand  side  of  (2.14)  can  be  rewritten  as 

OO 

y[n\  =  ^2  h{k)  [Kxmx]x[Kx  [n/Kx\  -  Kxmx\,  (2.15) 

mx=— oo 

where  |_  J  represents  the  floor  (integer  part)  of  its  argument,  and  the  term  Kx  -  [n/Kx\ , 
on  the  right-hand  side  of  the  equation,  is  needed  to  account  for  fact  that  x  is  not 
observed  at  every  time  n.  The  most  recent  observation  of  x ,  in  the  fundamental 
layer,  is  described  by  Kx[n/Kx\.  For  example,  if  Kx  is  3,  then  x  is  observed  only 
at  n  —  0,  3,  6, ....  Thus,  at  n  =  5,  for  instance,  the  most  recent  observation  of  x  is  at 
time  Kx[n/Kx\  =  3[5/3j  =  3  •  1  or  n  =  3.  Now  recall  that  a;  is  a  decimated  version 
of  x  that  satisfies 
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If  a  corresponding  decimated  version  of  the  filter  sequence  is  defined  as 

/r(fc)[  m  }  =  U»[KX  m], 

then  the  filter  equation  (2.15)  can  be  expressed  as 

OO 

y[n]  =  22  h{k)[mx]x[[n/Kx\  -mx\.  (2.16) 

mx=— oo 

Further,  observe  that  the  output  y[n]  has  values  valid  only  at  integer  multiples  of  Ky , 
i.e.,  n  =  rriyKy.  Thus  (2.16)  becomes 

OO 

y[Kymy\=  22  h{k)[mx\x[[Kymy/Kx\  -  mx\,  (2-17) 

mx=—o o 

Finally,  by  noting  that  y  is  a  decimated  version  of  y,  one  can  write 

(2.18) 

When  the  filtering  is  causal,  the  lower  limit  on  the  sum  in  (2.18)  can  be  changed  from 
mx  =  —  oo  to  mx  =  0  since  the  corresponding  impulse  response  terms  are  zero. 

Several  implementations  of  LPTV  filters  using  LTI  filters  are  possible. 
Two  common  representations  known  in  the  literature  are  shown  in  Fig.  2.13  and  Fig. 
2.14.  Figure  2.13  shows  an  LPTV  filter  realized  in  a  rotary  switch  form.  Given  a 
system  period  K,  the  rotary  switch  steps  through  the  sequence  0  <  k  <  K  —  1  in 
synchronization  with  the  rate  of  the  output  sequence  and  selecting  a  different  LTI 
filter  to  process  the  input.  While  this  form  is  typically  used  in  the  literature  to  apply 
when  the  filter  input  and  output  rates  are  the  same,  it  easily  applies  to  the  more 
general  case  represented  by  (2.18). 

Figure  2.14  shows  an  LPTV  filter  realized  in  filter  bank  form.  Here  the 
output  y[n]  is  equal  to  the  output  of  H^k\z)  where  n  =  k  (mod  K).  The  purpose 
of  decimation  followed  immediately  by  expansion  is  to  remove  any  non-zero  data  not 
associated  with  the  proper  phase  of  the  system  (for  a  detailed  explanation  of  this 
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Figure  2.13.  Rotary  Switch  Representation  of  an  LPTV  Filter 


Figure  2.14.  Filter  Bank  Representation  of  an  LPTV  Filter  (From  [Ref.  17]) 


filter  bank  see  [Ref.  17]).  The  Liter  bank  form  of  Fig.  2.14  is  restricted  to  situations 
where  the  input  and  output  are  sampled  at  the  same  rate. 

A  generalized  LPTV  Liter  bank  was  derived  that  allows  for  the  input 
x[mx]  to  be  observed  at  a  rate  Fx  and  the  output  y[my]  to  be  observed  at  a  rate  Fy. 
This  more  general  Liter  bank  is  depicted  in  Fig.  2.15.  This  Liter  bank  operates  at 
the  fundamental  rate  F,  thus  the  input  signal  x[mx]  needs  to  be  expanded  to  the 
fundamental  rate  by  Kx  and  the  output  needs  to  be  decimated  by  Ky.  In  addition 
the  delays  are  shifted  by  multiples  of  the  output  decimation  factor  Ky  and  there  are 
My  sets  of  Liters  H^ky\z),  where  My  is  the  number  of  samples  per  period  for  signal 
y  and  rny  =  ky  (mod  My).  If  Kx  and  Ky  are  coprime,  then  My  =  Kx]  but  this  is 
not  true  in  general.  Since  x[mx]  is  observed  at  a  rate  — —  of  the  clock  rate,  the  Liter 


coefficients  for  the  filters  H^ky\z)  are  of  the  form 

H(yky\z)  =  h0  +  h1z~Kx  +  h2z~2Kx  H - f  hP-iz~(p~1)Kx , 

where  P  is  the  filter  order. 

The  LPTV  filter  bank  of  Fig.  2.15  requires  a  clock  that  operates  at 
the  fundamental  rate  F.  If  two  clocks  are  available  that  operate  at  the  input  and 
output  rates,  then  the  LPTV  filter  bank  can  be  represented  as  shown  in  Fig.  2.16. 
In  this  form  z~ 1  =  z~Ky  and  z~Y  =  z~Kx .  The  decimation  and  expansion  factors  are 
replaced  with  My  since  the  right  hand  side  of  the  filter  is  operating  at  the  output  rate 
Fy.  In  addition,  the  Liters  can  be  described  by 

H^ky\zx)  =  h0  +  hiz-1  +  h2z~2  H - b  hP-iZ^p^ . 


C.  STATISTICAL  REPRESENTATION  OF  RANDOM 
SIGNALS 

In  order  to  adequately  develop  methods  and  models  for  processing  random 
signals,  it  is  important  to  be  able  to  describe  the  statistical  characteristics  of  the 
random  signals.  Since,  in  most  cases  the  first  and  second  moments  are  sufficient 
to  describe  the  necessary  statistical  characteristics,  this  section  is  limited  to  “second 
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clock  rate  =  F 


Figure  2.15.  Filter  Bank  Representation  of  an  LPTV  Filter:  Single  Clock 
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X 


Figure  2.16.  Filter  Bank  Representation  of  an  LPTV  Filter:  Dual  Clock 


moment”  properties.  It  will  be  seen  that  even  with  this  limitation,  there  are  significant 
issues  to  be  addressed  in  the  case  of  multirate  signals  and  systems. 

The  first  order  moment,  or  mean ,  of  a  random  process  is  defined  as  the  expec¬ 
tation  of  the  signal: 

lix[m\  =  £{x[m\}.  (2.19) 

The  mean  is,  in  general,  a  function  of  time,  although  for  wide  sense  stationary  (WSS) 
signals  it  is  constant.  For  cases  involving  multirate  signals  and  systems  the  mean 
is  (in  general)  a  periodic  function  of  time,  i.e.,  fix[m]  =  px[m  +  M\  where  M  is  the 
period. 

The  second  order  moments  of  a  random  process  are  defined  via  the  autocorre¬ 
lation  function.  Let  m  and  m'  represent  two  different  values  of  the  time  index.  Then 
the  traditional  autocorrelation  function  is  defined  as 

Rx[m,  m']  =  £{x[m\x*[m!]}  (2.20) 

where  V  denotes  the  complex  conjugate.  The  autocorrelation  function  is  used  for 
quantifying  second  order  statistical  linear  dependence  between  samples  of  the  signal 
at  different  points  in  time. 
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The  autocovariance  function  is  defined  as  the  autocorrelation  function  with 
the  mean  removed: 

C x[m,  m'}  =  £{(x[m]  -  nx[m])(x[m'}\  -  nx[m'])*}.  (2.21) 

Frequently,  the  autocovariance  is  more  useful  than  the  autocorrelation,  although  the 
two  functions  are  identical  when  the  mean  is  zero.  The  autocorrelation  and  the 
autocovariance  functions  are  related  by, 

Rx  [■ m ,  ml]  =  Cx  [■ m ,  m]  +  nx  [m]n*x [ml] .  (2.22) 

By  defining  the  term  lag  as  the  difference  between  the  two  time  instances, 
l  —  m  —  m',  the  correlation  and  covariance  functions  can  be  represented  also  in  terms 
of  the  time  index,  m,  and  the  lag,  /, 

Rx[m\  1}  =  £{x[m]x*[m  —  /]}.  (2.23) 

This  form  is  referred  to  as  the  “time-dependent”  autocorrelation  function,  or  the 
“time-lag”  autocorrelation  function.  An  analogous  definition  can  be  made  for  the 
autocovariance  function. 

Two  special  cases  are  important  for  multirate  systems.  First,  if  the  random 
process  is  stationary,  at  least  in  the  wide  sense,  then  the  mean  of  the  random  process 
is  constant  and  the  correlation  and  covariance  functions  are  dependent  on  the  lag 
only.  Thus  the  dependence  on  m  in  (2.23)  can  be  dropped  and  the  autocorrelation 
function  can  be  written  as: 


Rx[l]  =  £{x[m]x*[m  —  /]}. 

In  this  case  the  traditional  autocorrelation  function  Rx[m,  m'],  defined  in  (2.20)  de¬ 
pends  only  on  the  difference  l  =  rn  —  ml . 

A  second  case  of  importance  for  multirate  systems  is  the  cyclostationary  case. 
A  more  complete  discussion  of  cyclostationarity  can  be  found  in  [Ref.  31,  32],  however, 
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it  is  simply  observed  here  that  for  a  cyclostationary  process  the  time-lag  autocorre¬ 
lation  function  is  a  periodic  function  of  m,  i.e., 

Rx[m;  l ]  =  Rx[m  +  M;  l],  3 M  (2.24) 

where  M  is  the  period.  This  equation  will  be  taken  here  as  the  definition  of  a  cyclo¬ 
stationary  process.  Note  that  Rx[m\  l]  is  not  necessarily  periodic  in  /.  Although  for  a 
cyclostationary  random  process  the  traditional  autocorrelation  function  Rx  [m,  m']  is 
periodic  in  both  arguments,  with  the  same  period  M  (See  [Ref.  31,  32].) 

For  two  signals  sampled  at  the  same  rate  the  traditional  cross-correlation 
function  is  defined  as 

Rxy[m,m']  =  £{x[m\y*[m']}, 

where  m  and  m!  are  two  arbitrary  values  of  the  index.  The  cross-correlation  can  also 
be  represented  in  the  time-lag  form  as 

Rxy\m\  l ]  =  £{x[m\y*[m  -  /]}, 

where  l  —  m  —  m'.  Further  if  the  signals  are  jointly  wide-sense  stationary,  then  the 
cross-correlation  is  a  function  of  its  lag  term  only, 

Rxy[l]  =  £{x[m\y*[m  -  /]}. 

In  multirate  signal  processing  however,  one  must  be  able  to  find  the  cross-correlation 
of  signals  at  different  rates.  This  requires  more  careful  consideration.  The  discussion 
here  follows  the  development  in  [Ref.  35]. 

The  cross-correlation  function  for  signals  in  the  fundamental  layer  can  be  de¬ 
fined  in  two  forms: 

Rxyl^lxi  Tly\  £{x[7lx]z/  [/h/] }  (2.25) 

and 

Rxy[n;  v\  =  £{x[n]y*[n  -  v]}.  (2.26) 

The  latter  is  the  time-lag  form. 
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The  correlation  function  for  two  signals  at  different  rates  can  be  defined  in 
the  traditional  form.  Following  the  discussion  in  Section  II. B,  x[mx]  and  y[my ]  are 
represented  in  terms  of  signals  in  the  fundamental  layer.  This  representation  makes 
it  straightforward  to  write 

RXy[mx,my}  =  £{x[mx]y*[my]}  =  £{x[Kxmx]y*[Kymy\} 

Rxy  \R-xTnXi  -R y 'ixiy j .  (2.27) 

One  can  consider  the  cross-correlation  function  Rxy  [mx ,  my\  as  samples  of  the 
cross-correlation  Rxy  in  the  fundamental  layer,  as  shown  in  Fig.  2.17.  The  samples 
are  defined  on  a  rectangular  grid  or  lattice  [Ref.  59,  60]  of  points  separated  by  Kx 
samples  of  the  fundamental  layer  in  the  x  direction  and  K y  samples  in  the  y  direction 
(see  Fig.  2.17).  Correlation  for  values  of  nx  and  ny  not  on  this  lattice  are  not  defined. 


Figure  2.17.  Cross-correlation  Lattice. 

Developing  the  time-lag  form  of  correlation  for  two  signals  at  different  rates  is 
somewhat  more  complex.  Consider  the  following  definition,  which  can  be  written  by 
analogy: 

772  —  .  .  .  —  1  0  1 

Rxy[m;  l }  =  £{x[m]y*[m  -  l}}  =  £{x[Kxm\y*[Ky(m  -  /)]}, 

Z  =  . ..,-1,0,1,... 

(2.28) 
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With  v  as  the  time  delay  in  the  fundamental  layer,  the  cross-correlation  can  be  plotted 
in  the  (n,  u)  space  as  shown  in  Fig.  2.18.  The  indices  m  and  l  map  to  points  on  a 
non-rectangular  lattice.  Note  that  correlation  values  do  not  exist  for  all  possible 
combinations  of  n  and  v  because  samples  of  y  occur  only  when  its  argument  is  an 
integer  multiple  of  Ky. 


•  one  period 


v,l 


Figure  2.18.  Time-lag  Cross-correlation  Lattice. 


The  lattice  is  described  by  two  basis  vectors  Vi  and  V2  which  are  defined  in 
the  fundamental  layer  as 


Vl  = 


Kx 

Kt  -  K„ 


and 


v2  = 


0 

K„ 


The  indices  m  and  l  in  Rxy  [m;  /]  represent  a  point  on  the  lattice  that  is  reached  by 
taking  m  steps  in  the  Vi  direction  and  l  steps  in  the  v2  direction.  The  time-lag 
correlation  function  and  the  traditional  cross-correlation  function  are  related  as 


Rxy  ■>  ^  ]  Rxy\p~^t  rWl  ^]  •  (2.29) 

When  the  signals  x  and  y  are  jointly  wide  sense  stationary,  then  the  cross¬ 
correlation  as  defined  in  (2.28)  is  a  function  of  v  only  (and  not  n) .  Therefore,  the 
values  of  correlation  depend  only  on  the  distance  from  the  n  axis  in  Fig.  2.18  and 
since  the  sampling  pattern  is  periodic,  the  values  of  correlation  repeat  periodically. 
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Consider  the  following  example.  Suppose  the  cross-correlation  function  for 


signals  in  the  fundamental  layer  is  given  by 

RSy[nx,ny\=p-^~ny (\P\  <  1),  (2.30) 

where  nx  represents  the  time  instance  for  the  signal  x  in  the  fundamental  layer  and 
riy  represents  the  time  instance  of  y  in  the  fundamental  layer.  Note  that  the  symbol 
Rxy  is  used  to  distinguish  the  cross-correlation  of  the  signals  in  the  fundamental  layer 
from  Rxy,  the  cross-correlation  of  the  associated  decimated  signals. 

From  (2.27),  the  cross-correlation  of  the  decimated  signals  is 

Rxy  K,  my\  =  £{x[Kxmx\y*[Kymy]} 

Rxy  \KxWlx  f  KylTlyl 
ft—\Kxmx—Kymy\ 


Then  from  (2.29),  the  cross-correlation  can  be  represented  in  the  (m,  /)  space  as 


Rxy  [m- 1 }  =  Rxy  [m,  m  -  l] 

__  0—\Kxm—Ky(m—l)\ 
=  p  —  \  (Kx  —  Ky  )  m+Ky  1 1 


(2.31) 


Note  that  the  cross-correlation  of  the  multirate  signals  is  a  subset  of  the  values  of 
cross-correlation  of  the  signals  in  the  fundamental  layer. 

In  order  to  further  illustrate  the  relation  between  the  (mx,my)  and  (m,l) 
spaces,  the  following  example  is  provided.  In  this  example  numerical  values  are  used 
to  show  how  valid  data  points  in  the  (mx,my)  space  are  transformed  into  valid  data 
points  in  the  (m,  /)  space.  Assuming  Kx  =  2  and  Ky  =  3  and  using  the  definitions 
that  m  =  rnx  and  l  =  mx  —  my,  Table  2.1  shows  the  results  of  the  transformations 
and  Fig.  2.19  visually  depicts  how  the  points  map  from  one  space  to  the  other. 
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Table  2.1.  Index  transformation  of  the  cross-correlation 


(mx,my) 

(0.1) 

(1.0) 

(1.2) 

(2,1) 

(2,2) 

(  ■m;  l ) 

(0;  —1) 

(111) 
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(2;  l) 

(2;  0) 

Symbol 

0 

□ 

A 

0 

V 

•  *3 


EE 


A  ^ 
■©— 
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Figure  2.19.  Index  Transformation  of  the  Cross-correlation. 


D.  INPUT-OUTPUT  CROSS-CORRELATION  FOR  EX¬ 
PANSION,  DECIMATION  AND  FILTERING 

Having  defined  the  mean,  autocorrelation  and  cross-correlation  functions,  ex¬ 
pressions  can  be  derived  for  the  relationships  between  these  statistics  for  the  input 
and  output  for  some  common  operations.  In  particular,  the  operations  of  decimation, 
expansion  and  linear  filtering  are  considered.  In  the  following,  it  is  assumed  that 
the  autocorrelation  function  Rx[mx,m'x\  of  the  input  as  given  by  (2.20)  is  known. 
Equivalently,  the  time-lag  form  of  the  autocorrelation  Rx  [rnx  \  lx\  of  the  input  as  de¬ 
fined  by  (2.23)  is  known.  Further,  if  the  signal  is  wide-sense  stationary,  the  time-lag 
autocorrelation  function  can  be  reduced  to 

Rx[lx\  =  £{x[mx)x*[mx  -  lx\}.  (2.32) 
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1.  Decimation 

The  operation  of  decimation  is  depicted  in  Fig.  2.5  and  defined  as  (2.6).  For 
ease  of  reference,  the  mathematical  description  of  decimation  is  also  provided  below: 

y[my\  =  x[Lmy\  for  m  —  . . . ,  —  1,  0, 1, . . . . 

Given  two  arbitrary  time  instances,  rny  and  m'y,  the  output  autocorrelation 
for  the  decimator  is  given  by 

Ry  [my,  m'y]  =  £{y[my]y*[m'y\} 

=  £{x[Lmy\x*[Lm'y]} 

Ry  [my,  m’y]  =  Rx  [Lmy,  Lm'y\  (2.33) 

Thus  decimation  results  in  decimation  of  the  auto-correlation  function. 

By  defining  a  lag  term  which  is  the  difference  between  the  two  time  instances, 
l  =  rriy  —  'm'y  and  renaming,  m  =  myi  the  autocorrelation  function  can  be  represented 
in  the  time-lag  form  as 

Ry[m;  l }  =  £{y[rn]y*[rn  -  /]} 

=  £{x[Lm]x*[Lm  —  LI]} 

or 

Ry[m;  l \  =  Rx[Lm ;  LI]  (2.34) 

Further,  if  the  signal  is  wide-sense  stationary,  then  the  autocorrelation  function  is 
independent  of  m  and  can  be  written  as 

Ry[l\  —  Rx[Ll\  (2.35) 

The  cross-correlation  function  for  decimation  is  given  by 

RXy[mx,my]  =  £{x[mx]y*[my]} 

=  £{x[mx\x*[Lmy]} 
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or 


Rxyl'm'xi 'Wlyl  -Re  \RIxi  LlTly^  (2.36) 

Thus  the  original  cross-correlation  function  is  decimated  in  one  argument  only. 
Figure  2.20  visually  depicts  the  effect  of  decimation  on  the  cross-correlation  function. 
The  circled  locations  of  the  input  autocorrelation  Rx  are  the  values  preserved  in  the 
cross-correlation  of  the  decimator. 


Figure  2.20.  Relation  of  the  Input  Autocorrelation  to  the  Cross-correlation  of  the 
Decimator. 


Using  the  time-lag  definition  for  the  cross-correlation  of  two  signals  at  different 
sampling  rates,  the  cross-correlation  of  the  decimator  can  be  written  as 

Rxy\p^xi  ly\  £{x[771j;]?/  \tTLx  ly\~\ 

=  £{x[mx]x*[L(mx  -  ly )]} 

=  £{x[mx\x*[mx  -  ( Lly  -  (L  -  l)/^)]} 


or 

Rxylj^xi  ly\  R.r  ( h  l)^u]  (2.3T) 

The  expression  for  the  time-lag  cross-correlation  function  in  terms  of  the  auto¬ 
correlation  function  is  rather  cumbersome.  The  geometric  representation  in  the  (m,  u) 
is  more  clear.  The  effect  of  decimation  on  the  cross-correlation  in  the  (rn ;  l )  space  is 
visually  depicted  in  Fig.  2.21.  Again,  the  circled  values  of  the  input  autocorrelation 
are  preserved  in  the  cross-correlation. 
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Figure  2.21.  Relation  of  the  Input  Autocorrelation  to  the  Cross-correlation  of  the 
Decimator  in  the  (m;  l)  Space. 


Even  when  the  input  is  wide-sense  stationary,  the  cross-correlation  cannot  be 
described  by  its  lag  term  only,  since  the  lag  term  on  the  right-hand  side  of  (2.37) 
requires  knowledge  of  mx. 


2.  Expansion 

Again,  for  ease  of  reference,  the  mathematical  description  of  expansion  is  pro¬ 
vided  below: 


I  x[m/I]  m  div  / 
y[m\  =  l 

I  0  otherwise 

The  autocorrelation  of  the  output  signal  is  defined  as 


Ry[my,  rn'y\  =  £{y[my]y*[my}} 

£{x[my/ I]x*[m'y/ 1]}  my,m'y  div  / 

otherwise 


Therefore, 


Ry\my,m'y\  = 


Rx[my/I],m‘ '  /I]  myim'  div  / 


0 


otherwise 


where  my  and  m'  represent  two  arbitrary  values  of  the  time  index. 
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The  time-lag  form  of  the  autocorrelation  can  be  written  as 

Ry[my]ly]  =  £{y[my]y*[my  -  ly]} 

{£{x[my/I]x*[(my  -  ly)/I}}  my,  (my  -  ly)  div  I 
0  otherwise 

Therefore, 

Rx[my/I]-ly/I]  my,ly  div  / 

Rylmy >  y\  ~  ^ 

I  0  otherwise 

In  the  case  of  the  expander,  the  output  is  never  wide-sense  stationary  even 
if  the  input  is  WSS,  because  zeros  have  been  added.  Therefore,  the  autocorrelation 
of  the  expander  cannot  be  represented  by  a  lag  only,  and  a  time  reference  must  be 
stated  as  well.  If  the  input  is  WSS,  the  output  is  cyclostationary  however,  because 
the  autocorrelation  function  is  periodic  in  its  first  argument. 

The  cross-correlation  for  expansion  is  computed  as 


RXy[mx,my\  =  £{x[mx]y*[my}} 


{£{x[mx]x*[my/ /]}  rny  div  / 
0  otherwise 


Therefore, 


RXy  XflXl  Tfly 


I  Rx  [mx ,  my/I]  rriy  div  I 
I  0  otherwise 


Thus  the  original  cross-correlation  function  is  expanded  in  one  argument  only. 
Figure  2.22  visually  depicts  the  effect  of  expansion  on  the  cross-correlation  function. 
The  circled  locations  of  the  expander  cross-correlation  are  the  autocorrelation  values 
preserved  from  the  input.  All  other  values  are  zero. 
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Figure  2.22.  Relation  of  the  Input  Autocorrelation  to  the  Cross-correlation  of  the 
Expander. 


Using  the  time-lag  definition  for  the  cross-correlation  of  two  signals  at  different 
sampling  rates,  the  cross-correlation  of  the  expander  can  be  written  as 


Rxy[mx ;  ly\  =  £{x[mx]y*[mx  -  ly] } 

J  £{x[mx\x*[(mx  -  ly)/I}} 

lo 


(mx  -  ly)  div  / 

otherwise 


{£{x[mx]x*[mx  -  (ly  +  (/  -  1  )mx)/I)]}  (mx  -  ly)  div  I 
0  otherwise 


Therefore, 


\  Rx[mx,  (ly  +  (/  -  l)mx)/I)]}  mx-lyQ\ivI 
Rxy  \P^x 5  ly]  \  (2.38) 

I  0  otherwise 

The  effect  of  expansion  on  the  cross-correlation  in  the  (m,  l )  space  is  visually 
depicted  in  Fig.  2.23.  Again,  the  circled  values  of  the  cross-correlation  are  the 
values  preserved  from  the  input  autocorrelation.  The  uncircled  points  represent  valid 
locations  where  the  cross-correlation  is  defined,  but  these  values  are  zero. 

Similar  to  the  decimator,  when  the  input  of  the  expander  is  wide-sense  sta¬ 
tionary  the  output  cannot  be  represented  in  terms  of  its  lag  only.  The  first  reason 
is  that  the  lag  term  on  the  right-hand  side  of  (2.38)  requires  knowledge  of  mx.  In 
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Figure  2.23.  Relation  of  the  Input  Autocorrelation  to  the  Cross-correlation  of  the 
Expander  in  the  (m,l)  Space. 


addition,  the  expander  inserts  zeros,  therefore  the  output  is  not  WSS.  It  is  however, 
wide-sense  cyclostationary. 

3.  Linear  Time-invariant  Filters 

For  the  operation  of  linear  time-invariant  filtering  the  output  is  given  by 

OO 

y[m\  —  h[r\x[m  —  r]. 

r=— oo 

Given  two  arbitrary  time  instances,  m  and  m',  the  autocorrelation  of  the  filter  output 

is 


Ry[m,  m'}  =  £{y[m}y*[m'}} 


=  £  l  \  h[Ri]x[m  —  Ri]  J  I  h[R0\x[m'  —  R{ 

\Ri=— oo  /  \Ro=—oo 


0 


oo  oo 


CEE  h[Ri]h*[Ro]x[m  —  R\\x*{m'  —  R 

.  i?i=— oo  Ro=—oo 


0 


oo  oo 


-  E  E  h[Ri]h*[Ro]Rx[m  —  R±,m'  —  Ro 

Ri=—oo  Rq= — oo 

This  can  be  written  as  a  two-dimensional  convolution, 


Ry[m,  m]  =  R,x [m,  m]  *  h[m\  *  h*[m] 
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where  V  represents  convolution  with  respect  to  the  associated  variable. 

By  using  the  relation  Rx[m ;  /]  =  Rx[m,  m  —  l\,  where  l  =  m  —  m',  the  time-lag 
form  of  the  output  autocorrelation  becomes 

OO  OO 

Ry[m;  l]=  Y  h[Ri\h*[RQ}Rx[m  -  R^l  -  (Rx  -  R0)} 

i?l  =— OO  Rq  =—00 

For  a  stationary  process,  the  time-lag  correlation  function  is  independent  of  the  first 
argument.  Therefore,  the  result  for  stationary  processes  becomes 

OO  OO 

Ry[l]=  Y  Y  h[R-i\h*lRo\Rx[l  —  (i?i  —  i?0)]. 

R1=— oo  Rq=—oo 

This  can  be  put  in  the  well-known  form 


Ry[l]  =  Rx[l]*h[l}*h*[-l], 
where  V  represents  convolution. 

The  traditional  cross-correlation  function  between  the  input  and  the  output 
of  the  filter  is  given  by 


or 


Rxy[m,m']  =  £{x[m\y*[m']} 


—  £  <  x[m]  E  h[Ro]x[m  —  Ro 

l  \Rq=—oo 

{oo 

Y^  h*[Ro\x[m\x*[m'  —  Rq\ 

Rq=  oo 
oo 

=  Y  h*[Ro\Rx[m,  m'  -  R0\ 

Ro=—oo 


Rxy[m,  ml]  =  Rx[m,  m /]  *  h*[m!] 

Again  using  the  relation  Rxy[m;l\  =  Rxy[m,m  —  l],  where  l  —  m  —  m',  the 
time-lag  form  of  the  output  cross-correlation  function  is: 

OO 

Rxy  [m;  l]=  Y  h*  IR]  lm'J  +  ^o] 

Ro=—oo 
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or 

Rxy[m ;  /]  =  Rx[m;  l ]  *  h*[-l\. 

When  the  input  process  is  stationary,  this  reduces  to 


Rxy[l]  =  Rx[l]*h*[-l]. 


4.  Linear  Periodically  Time-varying  Filters 

When  calculating  the  correlation  functions  for  a  LPTV  filter,  care  must  be 
taken  to  account  for  the  difference  in  rates  between  the  input  and  output.  Recall 
that  the  general  equation  for  an  LPTV  filter  is  given  by  (2.18). 

The  traditional  cross-correlation  function  between  the  input  and  the  output 
of  the  LPTV  Liter  is  thus 


RXy[mx,my\  =  £{x[mx\y*[my\} 

=  £  lx[mx\  f  ^  h^k\R0]x[\rnyKy/Kx\  -  R0] 


\Ro=-oo 


£\  h*w[Ro}x[mx}x*[[myKy/Kx\  -  R0\ 


.  R0=-oo 


=  22  h*^[Ro]Rx[mx,[myKy/Kx\-R0]. 

Ro=—oo 

which  represents  a  convolution  for  h*(fc)  with  Rx  along  its  first  argument. 

The  autocorrelation  for  the  LPTV  Liter  is  given  by 

Ry[my,my]  =  £{y[my]y*[my}} 

=  £{(  12  h^lR^xHmyKy/K^-R,]]  (  J2  h^[R0]x[[m'yKy/Kx\  -  R0] 


oo  oo 


V  E  £  h^[Ri}h*^[R0\x[[myKy/Kx\  -  R1]x*[[m,yKy/Kx\  -  R0\ 

„  R\= — oo  Rq=—oo 


oo  oo 


V  £  h^[Ri}h^k,)[R0}Rx[lmyKy/Kx\  -  Ru  [m'yKy/Kx\  -  R0 ], 

R1=— oo  Rq  = — oo 
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where  my  =  k  (mod  Kx)  and  m'y  =  k!  (mod  Kx ).  This  equation  is  the  two- 
dimensional  convolution  along  indices  my  and  rriy. 

Because  of  the  floor  operation  used  to  transform  between  the  my  index  and  the 
mx  index,  the  time-lag  equivalents  for  the  autocorrelation  and  the  cross-correlation 
functions  become  rather  unwieldly. 

A  summary  of  correlation  relations  is  provided  in  Tables  2.2  through  2.4. 
A  summary  of  relations  in  the  frequency  domain  is  provided  in  Appendix  A.  These 
relations  were  not  used  in  this  dissertation,  so  are  not  presented  here,  but  are  provided 
in  the  appendix  for  reference. 

Table  2.2.  Summary  of  Decimation 
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Table  2.3.  Summary  of  Expansion 


AmA 


t  / 


II 

Rx[my/I ,  m'y/I]  my ,  m!y  div  / 

0  otherwise 

533 

«cS 

ii 

Rx[my/I]  ly/I\  my,  ly  div  I 
]  otherwise 

So 

H 

H 

II 

J  Rx[mx,my/I]  my  div  I 

1  0  otherwise 

S3 

H 

II 

Rx[mx]  {ly  +  {I  -  l)mx) / 1}  mx  —  ly  div  I 

0  otherwise 

Table  2.4.  Summary  of  Filtering 
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E.  MATRIX  REPRESENTATION 


In  order  to  develop  signal  processing  algorithms  in  the  time  domain,  it  is  fre¬ 
quently  easier  to  work  with  matrix  and  vector  representations  of  the  systems  and 
signals.  This  section  therefore  develops  tools  to  represent  the  basic  operations  of  dec¬ 
imation,  expansion  and  linear  filtering.  The  development  extends  the  ideas  described 
in  [Ref.  61]. 

1.  Decimation 

Consider  the  decimator  depicted  in  Fig.  2.5.  The  input  and  output  are  related 
by  y[m\  =  x[Lm\.  Define  the  two  vectors  of  input  and  output  samples 

r  i T 

x=  x[0]  x[l]  •••  x[LM  —  1] 

y=  2/[0]  y[  1]  •••  y[M  - 1]  • 

Then  the  decimation  operation  can  be  represented  as 

y  =  DLjMx, 

where  Dl,m  is  an  M  x  ML  matrix  of  l’s  and  0’s  used  to  extract  elements  of  x  as 
observations  to  be  placed  in  y.  The  matrix  D  l,m  will  be  called  a  decimation  matrix. 
For  example,  if  the  decimation  factor  is  L  =  3  and  the  size  of  the  desired  output 
vector  is  M  =  4,  then  the  associated  decimation  matrix  is 

100  000  000  000 
000  100  000  000 

D3,4  =  •  (2.39) 

000  000  100  000 

000  000  000  100 
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This  can  be  represented  in  a  more  compact  form  using  Kronecker  products.  In  this 
compact  form 

D  l,m  =  I  m  ®  d 

where  I m  is  the  M  x  M  identity  matrix  and  is  a  column  vector  consisting  of  a  one 
followed  by  L  —  1  zeros, 

dL  =  1  0  0  •  •  •  0  •  (2.40) 

Appendix  B  provides  a  short  discussion  of  Kronecker  products.  For  a  detailed  analysis 
of  Kronecker  products  and  matrix  calculus  see  [Ref.  62] ,  or  for  a  short  tutorial  paper 
see  [Ref.  63]. 

Using  this  definition  and  dropping  the  subscripts  for  ease  of  notation,  the  mean 
of  the  decimator  is 

m  y  =  E{  y}  =  £{Dx}  =  D£{x}  =  Dm^,  (2.41) 

and  the  output  autocorrelation  matrix  for  the  decimator  is 

R  y  =  £{yy*T}  =  £{Dx(Dx)*T}  =  D£{xx*T}D*T  =  DRSD*T.  (2.42) 

Likewise,  the  matrix  form  of  the  cross-correlation  is 

Rxy  =  £{xy*r}  =  ^{x(Dx)*r}  =  £{xx*T}D*T  =  RSD*T.  (2.43) 

For  two  separate  signals,  x\  and  X2,  that  are  decimated  as 

?/i  [m]  =  xi  [Lim]  2/2  [m]  =  x2[L2m], 
the  cross-correlation  matrix  is  given  by 

R.1,2  =  ^{yiy2T}  =  ^  L\,M\R-X\X2^*L2,M2  ■  (2.44) 

where  the  derivation  for  (2.44)  is  similar  to  that  for  (2.42)  above. 

The  covariance  matrix  is  defined  as 

C y  =  £{(y  -  m2/)(y  -  rrp,)*7’}  =  £{(Dx  -  Dmx)(Dx  -  Dmx)*T} 

=  D^{(x  -  ma.)(x  -  m2:)*T}D*:/ 
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or 

C  y  =  DCxD*r. 

The  following  relation  can  also  be  derived  easily, 


cv  =  =  D  (Rt  -  m„mJT)  D*T 


(2.45) 


2.  Expansion 

Now  consider  the  expander  defined  by  the  input-output  relation 

I  x[m/L\  m  div  L 
y[m]  =  l 

I  0  otherwise 

The  matrix  form  of  expansion  can  be  represented  as 

y  =  Ui^x, 

where  the  expansion  matrix  U  is  an  LM  x  M  matrix  of  l’s  and  0’s  used  to  insert 
L  —  1  zeros  between  samples  of  x  to  form  the  output  vector  y.  This  matrix  is  called 
an  expansion  matrix.  For  example,  if  the  expansion  factor  is  L  =  3  and  the  length  of 
the  input  vector  is  M  =  4,  then  the  associated  expansion  matrix  is 

10  0  0 
0  0  0  0 

0  0  0  0 

0  10  0 

0  0  0  0 

0  0  0  0 

U.3,4  = 

0  0  10 

0  0  0  0 

0  0  0  0 

0  0  0  1 

0  0  0  0 

0  0  0  0 
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The  expansion  matrix  can  be  represented  as 

U l,m  =  I m  ®  dL, 

where  I m  is  the  M  x  M  identity  matrix  and  d;  is  a  column  vector  defined  in  (2.40). 

Using  this  definition  and  following  the  procedure  used  for  decimation  in  the 
previous  section,  the  matrix  form  of  the  output  autocorrelation  for  the  expander  is 

R,  =  £{yy*r}  =  URSU*T. 

Likewise,  the  matrix  form  of  the  cross-correlation  is 

Rxy  =  £{xy*T}  =  RsU*t. 

The  matrix  results  are  of  little  use  by  themselves  because  there  are  large  blocks  of 
zeros  as  a  result  of  the  expansion  operation.  When  expansion  is  followed  by  linear 
filtering  however,  the  zeros  fill  in  so  a  more  meaningful  correlation  matrix  results. 
The  cross-correlation  of  two  interpolated  signals, 

I  Xi[m/Li]  m  div  L\ 

3/i  N  =  < 

I  0  otherwise 

{x-Am/L^  m  div  L2 
0  otherwise 

is  also  of  interest.  The  cross-correlation  matrix  for  these  two  signals  is 

^yiD2  =  £{yiY2  }  =  L1,M1R'XiX2^J  L2,M2' 

The  covariance  matrix  for  expansion  follows  a  similar  set  of  transformations. 
The  covariance  matrix  of  the  expander  is 

Cy  =  UCXU*T, 
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and  satisfies  the  relation 


C y  =  Ry  -  mym*yT  =  U  (Rx  -  nxrm*T)  U*T. 

A  direct  relationship  exists  between  the  decimation  matrix  and  the  expansion 
matrix.  The  expansion  matrix  is  the  transpose  of  the  decimation  matrix, 

U  l,m  =  D 

The  proof  of  this  result  follows  directly  from  the  definition 

U l,m  =  I m  <8)  dL  =  (lM  <S>  d^)  =  D£m, 

where  property  (3)  of  the  Kronecker  products  was  used  (as  described  in  Appendix 
B).  This  can  also  be  seen  in  the  examples  of  the  decimation  and  expansion  matrices, 
(2.39)  and  (2.46). 

Some  combinations  of  decimation  and  expansion  are  also  worth  mentioning. 
Expansion  followed  by  decimation  (using  the  same  factor)  is  represented  by  the  prod¬ 
uct  matrix  D iiMU LjM  =  D l,m m.  Since  this  pair  of  operations  results  in  no  change 
of  the  signal,  it  must  follow  that 

=  Im-  (2.47) 

On  the  other  hand,  decimation  followed  by  expansion  selects  every  Lth  term  of  a 
sequence  and  replaces  all  other  terms  with  zeros  (the  length  of  the  sequence  is  un¬ 
changed).  This  operation  is  represented  by  the  product 

=  I m  ®  did^.  (2.48) 
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This  product  matrix  is  an  orthogonal  projection  matrix.  The  geometric  interpretation 
of  this  operation  is  a  projection  of  the  vector  onto  a  subspace  defined  by  every  Lth 
component  of  the  original  vector. 


3.  Filters 

The  matrix  form  of  a  linear  filter  with  finite  support 

p+ 

y[m]  —  h[r\x[m  —  r]  (2.49) 

r=P_ 


can  be  represented  as 

y  =  htx, 

where,  the  term,  H,  represents  the  reversal  of  the  matrix  H.  which  is  found  by 
inverting  the  order  of  the  elements  along  both  the  columns  and  the  rows  (see  Appendix 
B.)1 

For  the  case  of  a  linear  time-invariant  filter,  the  matrix  FI  has  the  form 


h[P_\  0 

h[  1  +  P-\  h[P-\ 


0  0 
0  0 


H  = 


h[P+]  h[P+  -  1] 
0  h[P+\ 

0  0 

0  0 


0  0 
0  0 


h[P_]  0 

h[l+P-]  h[P- } 


0  0 
0  0 


h[P+]  h[P+  -  1] 
0  h[P+] 


1The  reversal  of  a  P  x  Q  matrix  A  with  elements  atj 


is  defined  as  a  matrix  A  with  elements 


aP—i,Q—j  • 
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where  P_  and  P+  are  the  values  of  the  lower  and  upper  limits,  respectively,  of  the 
filter  in  (2.49).  Using  this  definition,  the  matrix  form  of  the  output  autocorrelation 
of  the  filter  is 

R  y  =  £{yy*T}  =  HTR*H*. 

while,  the  matrix  form  of  the  cross-correlation  is 

R  xy  =  £{xy*T}  =  RaH* 

If  the  lower  limit  of  the  filter  matrix  P_  equals  zero,  then  the  filter  is  causal, 
since  the  filter  relics  only  on  past  observations  of  the  input  data.  On  the  other  hand, 
if  — P_  is  less  than  zero,  the  filter  uses  future  observations  of  the  input.  Therefore, 
the  filter  is  non-causal  and  cannot  be  implemented  in  real-time.  For  an  LPTV  filter 
there  are  up  to  K  possible  sets  of  filter  coefficients,  H^.  The  appropriate  set  of 
filter  coefficients  is  chosen  based  upon  the  system  phase,  as  discussed  earlier  in  this 
chapter. 

F.  SUMMARY 

In  this  chapter,  the  fundamental  building  blocks  and  necessary  indexing  schemes 
for  a  multirate  system  were  described.  In  addition,  methods  for  calculating  the  cor¬ 
relations  and  power  spectral  densities  for  these  building  blocks  were  presented.  Using 
these  building  blocks  and  relations,  it  is  possible  to  develop  the  multirate  Wiener- 
Hopf  equations  for  optimal  filtering  of  a  multirate  system.  This  is  the  topic  of  the 
next  chapter. 
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III.  LINEAR  ESTIMATION 


In  this  section,  the  methodology  for  optimal  filtering  with  multirate  obser¬ 
vations  is  developed.  The  basic  principles  of  linear  mean-square  estimation  apply 
to  multirate  signal  processing;  however,  due  to  the  periodic  nature  of  the  multirate 
system,  one  optimal  filtering  equation  is  not  sufficient  to  fully  describe  the  optimal 
filtering  problem.  By  exploiting  the  system  periodicity  of  the  multirate  system,  a 
set  of  optimal  filtering  equations  can  be  derived,  one  for  each  “phase”  of  the  filter. 
The  filter  output  (Wiener-Hopf  equations)  can  be  a  single  sequence  or  a  vector  of 
sequences  based  upon  the  input  sequences. 

A.  SYSTEM  EQUATIONS 

To  define  the  problem,  assume  a  multirate  observation  model  with  a  set  of  M 
wide-sense  stationary  observation  sequences,  Xi [mi]  ...  xm\p^m\-  These  sequences 
represent  M  observations  of  the  signal  s[n]  subjected  to  various  linear  degradations 
(e.g.,  additive  white  Gaussian  noise  (WGN)  or  linear  distortion),  shown  in  Fig.  3.1. 
These  observations  sequences  are  to  be  used  to  estimate  a  single  desired  sequence 
d[n],  as  shown  in  Fig.  3.2.  Let  F1,F2,.. .  ,Fm  represent  the  sampling  rates  of  the 
observation  sequences  and  let  K i ,  K2 , . . .  ,Km  be  the  associated  decimation  factors. 
It  will  be  assumed  that  the  estimate  d[n]  is  required  at  the  fundamental  rate  F,  so 
the  index  ‘n’  is  used  for  this  sequence.  The  desired  sequence  at  time  n  can  be  written 
as 

d[n]  =  d[K  ■  l  +  k] 

where  K  is  the  system  period,  l  =  \n/K\  and  n  =  k  (mod  K).  The  variable  k,  is  the 
sampling  phase,  as  defined  in  Chapter  II,  with  0  <  k  <  K  —  1. 
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Figure  3.1.  M-Channel  Multirate  Observation  Model 


Assume  that  a  causal  estimate  is  desired,  i.e. ,  the  estimate  is  calculated  using 
only  the  “present”  and  “past”  observations.  The  optimal  estimate  is  formed  by 
summing  the  output  of  linear  periodically  time-varying  (LPTV)  filters,  one  for  each 
observed  sequence  (see  Fig.  3.2).  Define  the  vector  of  filter  coefficients 


ht\Pi  - 1] 


1  <  i  <  M 


0  <  k  <  K  -  1 


where  P*  is  the  order  of  the  ith  filter  and  denote  the  vector  of  observations  from  the 
ith  channel  by 


Xi[mi]  Xi[rrii-  1] 


Xi[rrii  -  Pi  +  1] 


T 


where  r/i;  =  \_n/ Kt\ .  Here  m;  represents  the  most  recent  observation  time  for  ay  if 
the  causality  constraint  is  to  be  maintained.  Then  the  estimate  at  time  is  given 
by 

M  M 

dk[n ]  =  y^xf^r[n]  h\k ]  ^  h-fc)Tx-fc) [n];  where  n  =  k  (mod  K).  (3.1) 

1=1  2=1 

The  subscript  k  in  dk [n]  is  actually  redundant  since  n  =  k  (mod  K)\  however,  this 
notation  is  used  to  emphasize  that  the  statistical  properties  of  the  estimate  are  peri¬ 
odically  time  varying. 
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Figure  3.2.  Multirate  Optimal  Filter  (Direct  Form) 


The  form  of  the  estimation  given  by  (3.1)  and  depicted  in  Fig.  3.2  will  be 
called  the  “direct  form”  realization  since  the  input  sequences  are  directly  weighed  and 
summed.  In  the  next  section,  the  methods  for  finding  the  optimal  filter  parameters 
in  the  direct  form  and  the  corresponding  mean-square  error  variance  are  developed. 


B.  MULTIRATE  WIENER-HOPF  EQUATIONS 

Having  defined  the  estimate,  it  is  now  possible  to  define  the  error  as 

Ek[n]  =  d[n }  -  dk[n }  (3.2) 

and  to  find  the  optimal  set  of  filter  coefficients  that  minimize  £{|£fc[n]  |2}.  The  or¬ 
thogonality  principle  of  linear  mean-square  estimation  [Ref.  18,  58]  requires  that  the 
error  be  orthogonal  to  the  observation  vectors,  i.e. , 

£{xifc)M4M}  =  0;  ^  =  i,  2, ... ,  M-  (3.3) 

Substituting  (3.1)  and  (3.2)  into  (3.3)  and  taking  the  expectation  yields 

{/  M  \  *  M 

xfV]  yd[n]  -  JT[n]  hf]j  j  =  \f[n]  hf]*  =  0 
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where  the  following  terms  have  been  defined: 


hf  H  =  £{xfVK  W)  (3.4) 

R<‘)H  =  £{xf>[rl]xf*T[rl]}.  (3.5) 


This  can  be  rearranged  and  written  as 

M 

bf  =  i  =  l,2,...,A/ 

3=  1 


or  equivalently  as 


M 

^-ESTMhf.  (3.6) 

3=  1 

Then  using  the  Hermitian  symmetry  property  R^*  =  R^T,  (3.6)  can  be  represented 
in  matrix  form  as 


r  W 
Kn 

R  0)* 

R  W* 

hf)_ 

i 

CsT 1—1 

R  {k)T 

rv12 

R  (fc) 

rL22 

R  W* 

rv'2M 

= 

f(fc) 

1  d2 

f>(k)T 

_xvl  M 

R  (fc)T 

rV'2M 

R  (fe) 

f(fc) 

_ 

(3.7) 


The  associated  error  variance  according  to  the  orthogonality  principle  is  given  by 


crl  =  £{d[n]£*k[n\}  which,  upon  substitution  of  (3.1)  and  (3.2)  becomes 


a 


2 

k 


M 


R*( o)  -  E 

i= 1 


1  di 


(3.8) 


Observe  that  the  mean-square  error  is  periodically  time  varying,  i.e. ,  it  is  dependent 
on  the  phase  k  of  the  filter.  In  order  to  establish  a  single  figure  of  merit  for  the  system, 
an  average  error  is  defined, 

1 

/*'->  \  K 

b  =  n  b  <3-9) 

V  k= 0  / 
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This  geometrical  average  is  chosen,  since  the  error,  when  expressed  in  decibels  (dB) 
becomes  the  arithmetic  average  of  the  values  in  dB.  For  Gaussian  signals,  this 
relates  to  average  information  (entropy). 


C.  CALCULATION  OF  CORRELATION  TERMS 

The  correlation  terms  needed  in  (3.7)  and  (3.8)  can  most  easily  be  generated 
by  using  the  linear  algebra  concepts  described  in  Chapter  II.  For  any  time  n  and 
corresponding  index  k:  the  observation  sequence  x [n]  can  be  expressed  as 

xjfc)[n]  =  f>2)xi[n],  (3.10) 

where  x*  is  given  by 

r  i T 

Mn\  =  Xi[n]  Xi[n-  1]  •••  Xi[n  -  Pi  ■  lu  +  1\ 

and  D'^';1  is  an  appropriately-defined  decimation  matrix.  Note  that  x,  [n]  consists  of 
all  possible  values  of  x% [n]  (observed  or  unobserved)  and  x [n\  represents  just  the 
observed  values  used  by  the  filter. 

By  virtue  of  (3.10),  the  correlation  matrix  of  x, [n]  and  x^- [n\  is  given  by 

R«f  =  b£>  R,MT’  (3.11) 

where  Ry-  is  the  correlation  matrix  of  x,  [n]  and  x^  [n]  at  the  fundamental  layer 

R.tj  =  £{-ki[n\x*T[n]}. 

The  cross-correlation  vector  between  x,;  and  d  is 

b?  =  £>!%.  (3.i2) 
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where  r<&  is  the  cross-correlation  of  between  x*  and  d  at  the  fundamental  layer 

idi  =  £{Zi[n]d*[n\}. 

Note  that  the  dependency  on  the  system  phase  is  completely  determined  by  the 
decimation  matrices. 


D.  FILTER  COEFFICIENTS 


The  solution  to  the  multirate  Wiener-Hopf  equations  produces  the  parameters 
for  the  direct  form  realization  of  the  filter.  An  alternative  realization  can  be  developed 
that  leads  to  some  useful  insights  about  the  role  that  each  input  plays  in  forming  the 
overall  estimate  dfc  [n] .  This  alternative  realization  is  a  recursive  form  of  the  optimal 
filter  equations,  and  a  complete  derivation  of  these  equations  is  contained  in  Appendix 
C.  It  is  shown  there  that  the  filter  coefficients  for  the  direct  form  can  be  written  as: 

M 

if  ~  Gf'Tbf)  -  V  Ef>- 

j=i+ 1 


ppA)  _  Q(fc)*Tg(fc) 


,(*) 


(1  <i<  M  -  1) 


u(*0  _  -pAA1  U(k)  _  r(fc)*rr(fc) 
11 M  ~  1  dM  udM 


where  the  vector  b®  and  the  matrix  are  defined  as 


.(fc) 


(3.13) 

(3.14) 


1 

T— 1 

1 - 

'^5 

b(fc)  - 

udi  — 

B(fc)  = 
u 

=(fc) 

_  d(i— 1)_ 

L  h-bi J 
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and  the  matrices  G^  and  are  given  by 


G,W  =  { 


i  =  1 


R?° 

R(fc) 

rv12 

■p  (k) 

i o(k)*T 

B-12 

T>  (fc) 

rv2 

"R 

B(fc) 

ll 

R(fc)vr 

l/V-i) 

id  (k)*T 

1 

-l 


b ;j  1  <i<  m 


Ef ]  = 


R(n}  i  =  1 

Rife)  -  GW*TBlfc)  1  <  i  <  M 


The  hlter  coefficient,  h ^  in  (3.13),  can  then  be  re-expressed  as 

M 

h(k)  _  w(fc)  _  tt (fc)u(fc) 

i  i  h  j  ’ 


j=i+ 1 


(3.15) 


where 


and 


_ piW  1  ( cW)  (k)*7v  (fc) 

ni  ~  rdi  ~  Ddi 


JjW  _  gW  1 


With  these  definitions  it  can  be  seen  from  (3.15)  that  the  filter  h-/'’)  is  comprised  of 
a  modified  optimal  filter,  hfk\  and  cross  (or  prediction)  filters.  H-^,  which  allows 
the  optimal  filter  to  be  represented  in  the  form  shown  in  Fig.  3.3.  This  form  is 
referred  to  as  the  innovations  representation  because  each  branch  of  the  realization 
works  on  just  the  new  information  not  present  in  the  other  existing  channels.  It  can 
be  seen  that  after  applying  the  cross  filters  a  modified  set  of  observations  {x'±  •  •  •  x'M} 
is  obtained.  The  modified  observation  x\  is  the  residual  after  predicting  Xi  using  Xj, 
1  <  j  <  «  -  1  and  represents  the  new  information  brought  in  by  channel  i.  The 
primed  observations  are  mutually  orthogonal.  The  modified  optimal  filtering  terms, 
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h',  i  —  1,  2, . .  . ,  j,  represent  the  optimal  filter  for  estimating  d  if  Xj[rrij],  j  >  i  is  not 
used. 


(LPTV) 

Figure  3.3.  Multirate  Innovations  Representation 


By  separating  the  filter  into  cross  filter  terms  and  a  modified  optimal  filtering  term 
for  each  input  signal,  it  is  possible  to  determine  the  improvement  in  error  variance 
that  would  be  derived  from  adding  extra  signals,  without  having  to  recalculate  all  of 
the  filter  coefficients.  For  each  new  signal  added  only  the  prior  filtering  terms  and  the 
modified  optimal  filtering  terms  associated  with  the  new  signal  are  calculated.  All 
the  previous  prior  filtering  and  modified  optimal  filtering  terms  remain  unaffected. 

By  substituting  the  filter  coefficient  definition  of  (3.15)  into  the  error  variance 
of  (3.8)  the  error  variance  can  be  expressed  in  a  recursive  form  dependent  on  the 
number  of  signals  observed.  Let  al  ■  be  the  phase-periodic  error  variance  when  only 
j  signals  are  observed.  This  form  of  the  error  variance  will  be  called  the  innovations 
form  of  the  error  variance  equation  and  is  expressed  as 
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ak,M  ~  ak,M- 


M—l 

;(fe)T,  (k)*  _  ~(k)T 

cLM  11 M  1  clM 

1=1 

M—l  M 


u-(fe)*  \  u'  (k)* 

" niM  j  niM 


E  E 


( k )  T  ^ ( _ -pj(fc)*  ^  (k)*  _ 


LdM 


H 


i=  1  *i=i+l 
M—l  M  M 


hM  J 


1M 


M 


-EE  E  •••  E  ■ 

i=  1  il=i+l  12=*1  +  1  *M-l=*M-2+l 
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dM 
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(A?)* 
a  i 


x 


-H 


(fc)* 
'hi  2 


-H 


(fe)* 


-H 


*M- 1 


M 


(3.16) 


The  derivation  for  this  form  is  also  contained  in  Appendix  C.  Though  this  form  of 
the  error  variance  can  become  unwieldly  for  systems  with  a  large  number  of  observed 
signals,  this  equation  is  useful  to  show  the  improvement  in  performance  brought  about 
by  each  additional  channel. 


E.  PERFORMANCE  STUDY 

A  typical  scenario  for  multirate  estimation  is  based  on  the  signal  model  shown 
in  Fig.  3.4.  The  signal  of  interest  is  s[n] ,  therefore  d[n]  =  s[n].  The  desired  signal  is 
subject  to  independent  additive  noise  in  each  channel  before  decimation.  In  particu¬ 
lar,  this  study  considers  a  system  with  two  observation  sequences  Xi[mi\  and  .xq [m2] , 
with  decimation  factors  of  K\  =  2  and  K2  =  5. 

Two  different  types  of  signals  were  analyzed  for  comparison:  a  second  order 
AR  process  and  a  periodic  signal  comprised  of  two  sinusoids.  Both  signals  are  defined 
in  the  fundamental  layer  by  the  equations  shown  in  Table  3.1;  where  the  driving  term 
w\n]  in  the  AR  process  is  white  Gaussian  noise  (WGN)  with  mean  zero  and  unit 
variance. 

The  associated  autocorrelation  functions  for  each  of  the  signals  were  calculated  and 
are  provided  in  Table  3.2. 
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V\  [n] 


Figure  3.4.  Multirate  model  with  multiple  observation  sequences 
Table  3.1.  Signal  Examples 


Signal  Type 

Signal  Equation 

2nd  Order  AR 

s[n]  =  0.9s  [n  —  1]  —  0.8s[n  —  2]  +  Tu[n] 

2  Sinusoids 

s[n]  =  4cos[0.27m]  +  2cos[0.027rn] 

The  variances  of  the  two  noise  sequences  rji  and  7/2  were  chosen  such  that  the 
Signal  to  Noise  Ratio  (SNR)  assumed  the  values,  —6  dB,  —3  dB,  —1.7  dB,  0  dB,  1.7 
dB,  3  dB  and  6  dB.  The  corresponding  noise  variance  can  be  calculated  from 

2  _  Rs(°) 

aV  ]_q(SNR/10) 

where  Rs  is  the  autocorrelation  function  for  s[n].  A  block  diagram  of  the  optimal 
estimator  along  with  the  sampling  pattern  for  the  observations  is  shown  in  Fig.  3.5. 
Initial  simulations  were  conducted  with  filter  orders  of  P\  =  P2  =  5. 

Figure  3.6  shows  the  theoretical  error  variances  for  the  second  order  AR  model, 
computed  using  (3.8)  and  (3.9).  The  figure  compares  the  error  variances  for  using 
the  low-rate  signal  alone,  for  using  the  high-rate  signal  alone,  and  for  using  the  low- 
and  high-rate  signals  in  combination.  In  this  particular  experiment,  the  high-rate 
signal,  ,x'i ,  was  observed  in  a  0  dB  noise  environment  and  the  low-rate  signal,  X2  was 
observed  with  a  SNR  varying  over  the  range  —6  to  6  dB.  From  this  figure,  one  can 
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Table  3.2.  Signal  Autocorrelation  Functions 


Signal  Type 

Autocorrelation  Function 

2nd  order  AR 

Mi)  =  < 

3.7114  •  (0.8944)z  cos(1.0436Z  -  0.0646),  l  >  0 

-3.7114  ■  (1.1180)z  cos(1.0436Z  -  3.0770),  l  <  0 

2  sinusoids 

R s(l)  =  8cos(0.27tZ)  +  2cos(0.027tZ),  VZ 

*1  =M 


x2  =  [/«,] 


H 


■(*) 


dk  [n]  =  s  [«] 


K.  =2 


Channel  1 


Channel  2  ;o 


-o 


-cH 


P2 


K2=  5 


Figure  3.5.  Optimal  Estimator  Block  Diagram 


see  that  the  error  variance  for  the  low-rate  signal,  over  the  SNR  range  of  —6  dB  to 
+6  dB,  is  higher  than  the  error  variance  of  the  high  rate  signal  observed  at  0  dB. 
Individually,  the  high-rate  signal  performs  better  than  the  low  rate  signal.  When  the 
low-rate  and  high-rate  signals  are  processed  together,  however,  the  error  variance  is 
lower  than  for  either  signal  individually.  At  the  lowest  SNR  simulated  (—6  dB)  only 
a  modest  improvement  was  observed.  As  the  SNR  of  the  low  rate  signal  increases, 
however,  the  improvement  in  the  error  variance  approaches  2.5  dB. 

Figure  3.7  shows  the  theoretical  error  variances  for  the  sinusoidal  model. 
Again,  the  high-rate  signal,  was  observed  in  a  0  dB  noise  environment  and  the 
low-rate  signal,  X2  was  observed  with  a  SNR  varying  over  the  range  —6  dB  to  +6 
dB.  Again  in  this  case,  the  error  variance  for  the  low-rate  signal  at  all  SNRs  tested 
is  higher  than  the  error  variance  of  the  high-rate  signal  at  0  dB.  When  the  low-rate 
and  high-rate  signals  are  processed  together,  however,  the  error  variance  is  lower  than 
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2nd  Order  AR  Process:  High  Rate  SNR  =  0  dB 


Figure  3.6.  Error  Variance  vs  SNR  for  the  2nd  Order  AR  Process 

for  either  signal  individually.  At  the  lowest  SNR  simulated  (—6  dB)  only  a  modest 
improvement  was  observed.  As  the  SNR  of  the  low-rate  signal  increases,  however,  the 
improvement  in  the  error  variance  approaches  6  dB. 

Additional  simulations  were  conducted  with  filter  orders  Pi  and  P2  ranging 
from  5  to  30  for  an  environment  where  the  SNR  for  both  the  low  rate  and  high  rate 
signals  is  6  dB.  Selected  results  are  listed  in  Table  3.3  for  the  signal  with  two  sinusoids, 
showing  the  performance  that  results  when  using  both  the  high  rate  and  the  low  rate 
observations  together  or  using  only  one  set  of  observations  at  a  time.  From  Table  3.3, 
one  can  see  that  the  use  of  both  sets  of  observations  results  in  significant  improvement 
(3  dB  to  5  dB)  over  using  either  x\  or  X2  separately,  although  the  error  associated 
with  filtering  using  only  X2  is  large  compared  to  that  resulting  from  using  only  x\. 

The  effect  on  order  for  the  AR  process  could  not  be  studied.  Because  of  the  low 
order  of  the  AR  process,  varying  the  filter  orders  between  5  and  30  has  no  appreciable 
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Two  Sinusoidal  Input:  High  Rate  SNR  =  0  dB 


Figure  3.7.  Error  Variance  vs  SNR  for  the  Periodic  Signal 
effect  on  the  error  variance. 

F.  SUMMARY 

This  chapter  addressed  the  problem  of  optimal  filtering  of  multiple  related 
channels  of  data  observed  at  different  sampling  rates.  The  direct  form  of  the  Wiener- 
Hopf  equations  using  linear,  periodically  time-varying  filters  was  developed,  and  the 
error  variance  of  the  optimal  filter  was  observed  to  be  periodically  time-varying.  In 
addition  to  the  direct  form,  an  innovations  form  of  the  multirate  Wiener-Hopf  equa¬ 
tions  was  presented.  This  recursive  form  provides  insight  into  the  relative  change 
in  performance  when  additional  signals  are  added  or  removed  from  the  optimal  fil¬ 
ter.  Using  the  geometric  average  of  the  periodic  error  variances  as  a  performance 
measure,  it  was  demonstrated  that  optimal  Eltering  of  multiple  channels  can  provide 
improved  performance  over  optimal  filtering  using  one  channel,  even  if  the  secondary 
channel  when  used  alone  has  high  error  variances.  In  the  simulations  conducted,  the 
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Table  3.3.  Error  Variances  for  Periodic  Signals  with  2  sinusoids  (All  values  in  dB) 


Filter  Orders 

w2 

x i  [n]  and  X2[n] 

w2 

x  i  [n]  only 

w2 

X‘2  [n]  only 

P  =  5,Q  =  5 

-0.43 

3.27 

13.98 

P  =  30,  Q  =  5 

-13.54 

-12.26 

13.98 

P  =  5,Q  =  30 

-6.54 

3.27 

9.39 

P  =  30,  Q  =  30 

-17.32 

-12.26 

9.39 

largest  improvements  occured  when  the  signal  to  be  estimated  consisted  of  sinusoids 
in  noise.  When  the  signal  was  derived  from  a  second  order  AR  process,  the  use  of 
multiple  channels  and/or  higher  order  filters  resulted  in  only  a  small  improvement  in 
performance.  It  was  observed  that  a  single  channel  filter  of  order  5  was  sufficient  to 
estimate  the  second  order  AR  signal  in  noise.  These  examples  used  only  two  obser¬ 
vation  sequences,  but  the  derivations  and  experiments  presented  here  can  be  applied 
to  an  arbitrary  number  of  channels. 
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IV.  LINEAR  PREDICTION 


The  previous  chapter  discussed  estimation  of  data  from  noisy  signals  using 
present  observations  and  a  finite  number  of  past  observations.  Applications  that 
require  the  ability  to  estimate  the  present  value  of  the  observations  themselves  using 
only  past  observations  (causal  filtering),  however,  exist.  A  few  of  the  most  common 
applications  are  focused  on  speech  and  image  coding,  target  tracking  and  target 
classification.  The  process  of  estimating  a  signal  using  only  its  past  observations  is 
known  as  linear  prediction.  Linear  prediction  is  at  the  core  of  all  linear  estimation 
problems,  such  as  the  more  general  Wiener  filter.  Furthermore,  specific  insights  that 
arise  from  linear  prediction,  such  as  the  lattice  forms  of  the  filter,  can  be  carried  over 
to  the  more  general  problems. 

A.  SINGLE-CHANNEL  AND  MULTICHANNEL  RE¬ 
VIEW 

Linear  prediction  has  been  well  researched  for  the  single-channel  and  multi¬ 
channel  cases.  At  its  most  intrinsic  form,  the  goal  is  to  estimate  the  present  value  of 
a  data  sequence  from  a  finite  collection  of  past  data.  For  the  single-channel  case,  the 
prediction  equation  can  be  written 

x[n]  =  —a\x[n  —  1]  —  a*2x[n  —  2]  —  •  •  ■  —  a*Px[n  —  P } 

p 

=  ~^2a*x[n-i] 

i= 1 

where  x[n]  is  the  predicted  value  of  x[n]  using  x[n  —  1]  through  x[n  —  P]  and  a\ 
through  aP  are  the  prediction  filter  coefficients  (written  as  conjugate  negative  values 
for  later  convenience).  The  value  P  is  referred  to  as  the  prediction  filter  order. 
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For  the  multichannel  case,  the  prediction  equation  generalizes  to 
x[n]  =  —A*7  x[n  —  1]  —  A^1  x[n  —  2]  —  •  •  •  —  ApTx[n  —  P } 

p 

=  ~  AiT^ln  ~  *1  (4-1) 

i=l 

where  the  vector  x  represents  the  data  observed  in  each  of  M  channels  and  A]  through 
A  P  are  the  prediction  filter  coefficient  matrices.  This  more  general  case  will  now  be 
discussed. 

The  goal  of  the  linear  prediction  problem  is  to  find  an  optimal  solution  for 
the  coefficients  A,  in  (4.1),  such  that  the  expected  value  of  the  squared  norm  of  the 
prediction  error  e  [n]  =  x[n]  — x[n]  is  minimized.  The  solution  to  the  linear  prediction 
problem  results  in  the  well  known  (multichannel)  Normal  equations 

R[0]  R[l]  R  [P]  I  E 

R[-l]  R[0]  R[P-1]  A,  =  0 

R|-P]  R[— P  +  1]  R[0]  A  P  0 

where  the  multichannel  correlation  function  is  defined  as  R[i]  =  £{x.[n]x.*T[n  —  ?’]} 
and  E  is  the  prediction  error  covariance  matrix  E  =  £{e  [n]s*T[n]}. 

If  the  order  of  the  prediction  filter  is  allowed  to  grow  at  each  successive  ob¬ 
servation,  in  order  to  include  all  past  observations  at  each  step,  then  the  calculated 
prediction  errors  £ [n]  will  be  orthogonal.  If  a  prediction  filter  of  sufficiently  high 
order  is  used,  then  the  calculated  prediction  errors  will  be  approximately  orthogonal 
(£{s[i]£*T[j]  =  0}  for  i  ^  j).  Thus  the  prediction  error  filter  can  be  thought  of  as 
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a  causal  whitening  filter.  This  property  is  important  in  autoregressive  (AR)  modeling, 
moving  average  (MA)  modeling  and  autoregressive  moving  average  (ARMA)  modeling 
(see  [Ref.  18,  58,  19].) 


B.  MULTIRATE  LINEAR  PREDICTION  THEORY 


The  method  for  extending  linear  prediction  to  a  multirate  system  will  now  be 
considered.  Recall  from  Chapter  II,  that  multirate  systems  are  periodic  in  nature  and 
that  one  can  determine  an  associated  system  period  K .  This  system  period  can  be 
used  to  partition  the  data  into  blocks  with  a  common  sampling  structure  as  illustrated 
in  Fig.  4.1 

Channel  1 

Channel  2 

time 


o- 

/ 

t 

H-O— 

— O 

4c~ 

- o — • — 

Figure  4.1.  Multirate  System  Block  Structure 


In  this  example,  Channel  1  is  observed  at  half  the  fundamental  rate,  i.e. ,  the 
decimation  factor  is  2,  while  for  the  Channel  2  the  decimation  factor  is  3.  The 
system  period,  as  defined  in  Chapter  II,  for  this  example  is  K  =  6.  This  has  been 
highlighted  in  Fig.  4.1  by  the  dashed  boxes  representing  blocks  of  data.  When  viewed 
at  the  fundamental  rate,  the  statistics  of  the  multirate  system  are  cyclo-stationary. 
However,  when  the  multirate  system  is  viewed  in  blocks,  the  block  statistics  are 
stationary. 

In  the  specific  linear  prediction  problem  that  is  considered  here,  each  new  data 
point  (in  any  channel)  is  predicted  as  it  occurs  in  time.  The  prediction  is  based  on 
a  finite  number  Q  of  previous  full  blocks  of  data  as  well  as  all  of  the  data  within  the 
current  block  that  occurs  earlier  than  the  point  being  predicted.  Specifically,  let  m 
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represent  a  time  index  for  the  blocks  and  i  =  1,2, ...  ,P  be  an  index  for  the  data 
occurring  within  a  block,  also  ordered  in  time.  For  the  moment,  assume  that  data 
points  from  different  channels  do  not  occur  at  exactly  the  same  time.  The  prediction 
equation  for  the  ( i  +  l)st  data  point  is  then  written  as 


Xi[m]  =  —  a^Xi-ilm]  —  ■  ■  ■  —  a*i_1xi[m]  —  a*Jx[m  —  1]  —  ...  —  a  *qx[m  —  Q\ 

i—  1  Q 

=  -  atpxi-P N  “  a*J x[m-q].  (4.3) 

P=  1  9=1 


The  variable  Xi[m ]  represents  the  ith  observed  data  point  for  the  mth  system  block, 
and  x[m]  is  a  vector  of  all  the  observed  data  associated  with  the  rnth  system  block. 
The  elements  of  x[/n]  are 

xp[m] 


x[m] 


xP-i[m\ 


(4.4) 


xi[m } 

This  method  of  ordering  the  block  data  points  is  for  indexing  purposes  only;  it  does 
not  account  for  data  being  observed  at  the  same  time.  The  data  is  ordered  according 
to  oldest  data  first  with  channel  1  being  considered  the  oldest  channel  for  indexing 
purposes.  Figure  4.2  shows  these  variables  for  the  system  depicted  in  Fig.  4.1. 


Channel  1 


Channel  2 


t[m  —  l] 


x[»j] 


x2[m-\]  x4[m-\]  x2[m]  x4  [m] 


Figure  4.2.  Multirate  System  Block  Variables 


The  prediction  error  for  Xi[m]  is  given  by 

£i[m]  =  Xi[m\  —  Xi[m ] 
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or 


i— 1 


Q 


(4.5) 


Ei[m]  =  Xi[m]  +  V  a)iPXi_p[m]  +  2j  a*Jx[m  -  5] 

P=  1  9=1 

By  applying  the  principle  of  orthogonality  [Ref.  18],  one  can  determine  the 
prediction  error  variance  for  each  predicted  value  and  the  set  of  equations  that  the 
prediction  coefficients  must  satisfy.  The  prediction  error  variance  satisfies 


i—  1 


Q 


o]  =  £{xi  [m\e* [m] }  =  RXi  [0]  +  ^  aitPRXiXi_p  [0]  +  aI,grxiX [q]  (4.6) 

P= 1  9=1 

where  rXiX[l ]  is  the  vector  of  cross-correlation  terms 


:[*]  = 


RxiXp  [^] 
RxiXp- 1  [l] 


RxiXl  [l] 

Setting  the  observations  orthogonal  to  the  prediction  error  produces  the  equations 


0  =  £{xj[m  —  Z]£*[m]} 

p  Q 

=  RxjXi[— l]  +  y  p  ai,pRxjXj-p  [~^]  +  y  ^  ai,qrXjTclg  ~  l\- 
P=  1  9=1 


(4.7) 


Combining  these  sets  of  equations  for  all  observations  within  a  data  block,  and 
using  the  shorthand  notation 


R;  =  RX[Z]  =  £  [: x[m]x*T[m  —  /]]  , 
produces  the  matrix  form  of  the  multirate  Normal  equations 


(4.8) 


Ro  Ri 

R-i  Ro 

R-q  Ri-q 


Rq 

Rq-i 

Ro 


Ao 

E 

Ax 

= 

0 

_A  Q_ 

0 

(4.9) 
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where 


0  •••0 

1  •••0 

.-1,1  ■■■  o  (4.10) 

|_«P,P- 1  ttp— 1,P— 2  •  •  •  1 

Aj  elp,*  3-P— l,i  '  '  '  &l ^  •  •  •  i  Qi  (4,11) 

and  the  error  covariance  matrix  is 

CTp  x 

0 

0  0 

The  matrix  elements  represented  by  an  ‘  x  ’  do  not  affect  the  linear  prediction  problem, 
so  their  values  do  not  have  to  be  calculated. 

So  far,  it  has  been  assumed  that  data  points  from  different  channels  do  not 
occur  at  the  same  time.  If  points  do  occur  at  the  same  time  then  they  must  be  jointly 
predicted.  For  instance,  for  the  system  depicted  in  Figure  4.2,  data  points  occur  in 
both  channels  at  the  first  observation  time  of  the  data  block.  The  prediction  equation 
for  X\ [m]  is  of  the  form 

Q 

xi[m]  =  —  a i^x[m  —  q\. 

9=1 


1 

«p,  i 

A-0  —  (X  P2 
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This  equation  does  not  change.  The  linear  prediction  equation  of  (4.3)  for  X2  [m], 
however,  would  change  from 


x2[m]  =  -a*2^xi[m\  -  ^  a^x[m  -  q] 

9=1 


x2[m]  =  -  a2^x[™  -  ?]• 

9=1 


In  other  words,  xi[m]  and  x-i  [rn]  are  predicted  simultaneously.  The  prediction  coef¬ 
ficient  matrix  and  prediction  error  matrix  of  (4.10)  and  (4.12)  for  the  5x5  system 
with  no  joint  observations  are 


An  — 


1 

0 

0 

0 

0)5,1 

1 

0 

0 

Ot  5,2 

0)4,1 

1 

0 

«  5,3 

CD  4,2 

OL  3,1 

1 

0)5,4 

CD  4,3 

0)3,2 

Q)  2,1 

^5 

X 

X 

X 

X 

0 

crl 

X 

X 

X 

0 

0 

a3  X 

X 

0 

0 

0  a22 

X 

0 

0 

0 

0 

(4.13) 


(4,14) 


These  matrices  would  change  to 


An  — 


1 

0 

0 

0 

0 

Q)  5,1 

1 

0 

0 

0 

0)5,2 

0)4,1 

1 

0 

0 

0)5,3 

0)4,2 

Q)  3,1 

1 

0 

0)5,4 

0)4,3 

0)3,2 

0 

1 

(4.15) 
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and 


of  x  x  x  x 

0  of  x  x  x 

E  =  0  0  of  x  x  (4.16) 

0  0  0  of  of, 

0  0  0  of2  of 

because  of  the  joint  observations.  Notice  that  the  lower  right  block  in  A0  has  been 
changed  to  the  2x2  identity  matrix  while  the  lower  right  block  E  is  changed  to  a  2  x  2 
error  covariance  matrix  which  characterizes  the  joint  prediction  of  x,  and  x2.  The 
equations  for  solving  the  Normal  equations  must  be  slightly  modified  to  account  for 
any  joint  observations.  In  particular,  the  covariance  elements  of  the  joint  observations 
can  be  calculated  from 

^  Q 

Q ij  =  =  RxiXj  [0]  +  'y  ]  Otj,pRxjXi-v  [0]  +  y  ^  a-i,qrXix[<l\  >  (4.17) 

P  9=1 

where  i  and  j  are  jointly  observed,  the  summation  over  p  does  not  include  any  joint 
observations  and  of  =  erf.  In  addition,  the  orthogonality  equations  become 

0  =  £{xj[m  —  /]£*[m]} 

Q 

—  RxjXi[— l]  +  y  '  ai,pRxjXj-p [~^]  T  y  '  a.iQrXjX\q  —  /].  (4-18) 

P  9=1 

The  explicit  form  of  the  multirate  linear  prediction  problem  can  most  easily 
be  seen  by  using  an  example.  Figure  4.3  shows  a  two-channel  multirate  system  with 
decimation  factors  of  K1  =  1  and  K2  =  2.  Channel  1  has  two  samples  per  period  while 
channel  2  has  one  sample  per  period.  In  addition,  observations  from  both  channels 
occur  simultaneously  at  every  other  system  clock  interval.  Thus  joint  prediction 
of  the  channels  must  occur  at  those  times.  The  block  variables  xt  [rri]  used  in  this 
example  have  been  labeled  in  Fig.  4.3.  Using  a  filter  order  of  P  =  1  and  assuming 
that  the  block  of  observations  at  m  —  0  are  available,  the  linear  prediction  equations 
for  predicting  x% [m]  can  be  determined.  For  Xi[l]  the  equations  would  be 
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n 


Channel  1 

Channel  2 


xi  [2] 

— o 

X2  [0]  x2  [l]  x2  [2] 


m =0  m = 1 


o — •— 

*i[0]  -p[°] 

n  . 

• — o— 

y  [i]  X3  [i] 

kJ  • 

kJ  * 

Figure  4.3.  Illustration  of  Decimation  for  a  Linear  Prediction  System. 


xi[l]  =  -a{tlx3 [0]  -  a*h2x 2[0]  -  ai  3xi[0] 

The  error  would  then  be  £i[l]  =  aq[l]  —  xi[l]  or 

£i[l]  =  %i  [1]  +  [0]  +  a*lj2x2[  0]  +  ai)3xi[0] 

Applying  (4.17)  with  i  —  j  —  1,  the  prediction  error  variance  erf  would  be 

erf  —  Rx  1  [0]  +  CH,lRXlx3  [1]  +  ®1,2-Rxix2  [1]  +  ®1,3-Rxi  [1]  (4.19) 

and  the  prediction  error  covariance  for  observation  a^fl]  and  xi[l]  is 

a2i  =  Rx2x i [0]  +  a\y\RX2X3\\-]  +  ai,2-Rx2[l]  +  ®i,3-Rx2xi  [1]  (4.20) 

Applying  (4.18),  the  orthogonality  equations  for  the  previous  observations  are 

-Rx3xi[— 1]+  ®i,i-Rx3  [0]+ni,2-Rx3x2  [0]+ni,3-Rx3xi  [0]  =  0  (4-21) 

Rx2x\  [  1]  Toiii-Rx2x3  [0]  T  oi,2-Rz2[0]Tai,3-Rx2iri  [0]  =  0  (4.22) 

RXl  [- 1]  +oi,i-Rxix3  [0]  +oi,2-Rxix2  [o]  +  o>i}3RXl  [0]  =  0  (4.23) 

In  a  similar  fashion,  the  linear  prediction  equations  for  x2\\-}  can  be  written. 
The  prediction  equation,  using  (4.3),  is 

X2[l]  =  -«2,lXl[0]  -  a22X2[0]  -  a2  3Xi  [0] 
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with  an  error  equation,  from  (4.5),  of 


e2[l]  =  X2[l]  +  02,1X3  [0]  +  02,2X2(0]  +  02,3X1  [0] 

The  prediction  error  variance  and  prediction  error  covariance  equations,  derived  from 


(4.17),  are 

<x2  —  -Rx2[0]  +  02,l7?z2x3  [1]  +  02,2-Rz2[1]  +  0'2,3-Rx2o;i  [1]  (4.24) 

0,2  Rx\X2  [0]  T  02, 1^,3  [1]  T  Oj2  2RXix2  [1]  T  Oj2  ^RXi  [l]  (4.25) 

The  orthogonality  equations  for  the  previous  observations,  using  (4.18),  are 

Rx3x2  [— 1]+  a2,iRx3[0}+a2,2RX3X2[0]+a2}3RX3Xl[0]  =  0  (4.26) 

Rx2  [—  1]+02,1-Rx2x3  [0]+  0'2,2-Rx2  [0]+O:2,3-Rz2xi  [0]  =  0  (4.27) 

Rx3X2  [—  1]+02,1-Rxix3  [0]+O2,2-R*ia:2  [0]+  Oj2^RXi  [0]  =  0  (4.28) 

Finally,  the  equations  for  X3[l]  can  be  written  as 

x3[l]  =  -«3,iX2[l]  -  a3,2xi[l]  -  03,1X3  [0]  -  03,2X2  [0]  -  03,3X1  [0] 


£3[1]  =  X3[l]  +  «3,iX2  [1]  +  «3,2^1  [1]  +  03,1X3(0]  +  03,2X2(0]  +  03,3X1(0] 

and 

O3  Rx3  [0]  T  C^3,lRX3X2  [0]  T  (^3,2Rx3xi  [0]  T  &3,lRx3  [1]  4“  &3,2Rx3X2  [1]  T  ®j3,3Rx3xi  [l]  (4.29) 

There  are  no  prediction  covariance  equations  since  there  are  no  observations  occurring 
at  the  same  time  as  X3[l].  In  this  case  however,  the  number  of  orthogonality  equa¬ 
tions  increases  to  five,  since  the  error  £3[1]  must  be  orthogonal  to  all  five  previous 
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observations  (see  Fig.  4.3)  These  orthogonality  equations  are: 


RX2X3  [0]T  ^3,1-^o:2  [0]  ^3,2-^a:2^i  [^]  4_n3?i-R^2^3  [l]  T  ^3,2^a?2  [l]  Tn3,3-RtC2£i  [^]  0 

(4.30) 

^x±X3  [0] T  ^3,l-^xia;2  [0] 4"  ^3,2-Rtci  [0] ~\~CL^iRXiX3  [l] Tn3,2-Ra?i:r2  [l] T  ^3,3-ftci  [1]  0 

(4.31) 

Rx3  [— 1]+«3,1-R  X3X2  [— l]+a3,2-R  £3^1  [  — 1]+  03,1^0:3  [0]+Ct3,2-R  X3X2  [0]  ^U'i.'.iR-XzXl  [0]  0 

(4.32) 

Rx2X3  [  l]-!-  0!3^RX2\  1] H“C^3, 2RX2X1  [  1] ~ l- CI3, \Rx2X3  [0] 4”  ®3, 2Rx2  [0] 4“®3, ‘iRx2X\  [0]  0 

(4.33) 

Rxix3  [  1]  4~  (%3,lRxiX2  [  1]  4"  CY^  ^Rxi  [  1]  4“  ®3,l4?xia;3  [0]  4"  0‘3,2Rx\X2  [0]  4~  ^3,3Rxi  [0]  0 

(4.34) 

By  combining  (4.19)  through  (4.34),  the  matrix  form  of  the  multichannel  Nor¬ 
mal  equations  (4.9)  can  be  produced  and  written  as 


where 

Ro  Ri 

R  =  .  (4.35) 

R-i  :  Ro 
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or,  more  explicitly, 


RX3  [0] 

RX3X2  [0] 

Rx  3x\  [0] 

:  RxA1] 

RX3X2  [1] 

Rx3 Xl  [1] 

RX2X3  [0] 

RX2  [0] 

Rx 2X\  [0] 

•  RX2X3  [1] 

Rx2[  1] 

Rx  2X1  [l] 

RxiX3  [0] 

RX \X2  [0] 

flxJO] 

RX\X3  [1] 

RX\X2  [1] 

RxA1] 

RX3X2  [~  1] 

RX3X1  [—  1] 

:  RX3  [0] 

RX3X2  [0] 

Rx 3xi  [0] 

RX2X3  [  1] 

Rx2[-  1] 

RX2X1  [  1] 

RX2X3  [0] 

Rx2[  0] 

Rx  2X1  [0] 

_RXlX3  [  1] 

RXiX2  [  1] 

RxA- 1] 

Rx \x3  [o] 

Rx \X2  [0] 

74i[0]  _ 

1 

0 

0 

«3,1 

1 

0 

Ao 

Oi  3,2 

0 

1 

Ax 

^3,1 

02,1 

Ol,l 

®3,2 

02,2 

Ol,2 

03,3 

02,3 

Ol,3 

and  the  error  matrix  is 


X 

X 

0 

02 

E 

0 

*^12 

^i2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

(4.36) 


(4.37) 
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C.  AN  EFFICIENT  SOLUTION  TO  THE  MULTIRATE 
NORMAL  EQUATIONS 


In  the  previous  section,  the  Normal  equations  for  the  multirate  linear  predic¬ 
tion  problem  were  developed  element  by  element.  Although  these  equations  can  be 
solved  directly,  it  is  desirable  to  develop  a  method  for  calculating  the  linear  predic¬ 
tion  coefficients  and  error  variances  more  efficiently.  Efficient  recursive  methods  for 
solving  these  equations  have  been  developed  for  both  the  single-channel  and  mul¬ 
tichannel  cases.  In  the  single-channel  case  the  method  is  called  the  Levinson  (or 
Levinson-Durbin)  algorithm  while  in  the  multichannel  case  it  is  known  commonly  as 
the  multichannel  Levinson  algorithm,  or  in  some  literature  it  is  called  the  Levinson- 
Wiggins-Robinson  (LWR)  algorithm  [Ref.  64],  In  the  following,  a  similar  procedure 
is  developed  for  the  multirate  case. 

The  equations  to  be  solved  are 


Ro 

Ri  •• 

Rq 

Ao 

E 

R  i 

Ro 

•  Rq-1 

Ax 

= 

0 

R-q 

Ri-q  • 

Ro 

_Aq_ 

0 

These  are  the  multirate  Normal  equations  of  (4.9),  repeated  here  for  convenience. 

To  begin  the  solution  to  (4.38),  consider  the  problem  of  predicting  an  entire 
block  of  data  at  once.  This  is  identical  to  the  multichannel  problem.  The  Normal 
equations  for  this  problem  are  of  the  form  (4.2)  and  can  be  written  as 
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(4.39) 
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Observe  that  the  correlation  matrix  in  (4.39)  is  identical  to  that  in  (4.38)  although 
the  variables  A)  and  E'  are  not.  If  (4.39)  is  right  multiplied  with  Ao  however,  and 
(4.38)  compared  to  (4.39),  it  can  be  seen  that  the  following  equalities  must  hold: 

A,  =  A;Ao,  i  =  l  ,...,Q  (4.40) 

E  =  E'Ao.  (4-41) 

Therefore,  if  solutions  to  (4.39)  and  the  matrix  A0  can  be  found,  then  the  desired 
multichannel  parameters  A,:  (for  i  —  1,  2, . . , ,  Q)  and  E  can  be  found  from  (4.40)  and 
(4.41). 

A  computationally  efficient  solution  to  (4.39)  can  be  found  using  the  LWR 
algorithm  (see  Appendix  D).  This  solution  produces  the  matrices  A),  where  i  = 
1,2, ...  ,Q  and  E\  To  find  Ao  from  (4.41)  the  forms  of  E  and  Ao  must  be  considered. 
In  particular,  E  is  upper  triangular,  e.g.,  (4.14),  when  there  are  no  simultaneous 
observations,  and  block  upper  triangular,  e.g.,  (4.16),  when  some  observations  occur 
simultaneously.  Also,  A0  is  lower  triangular  with  unit  diagonal,  e.g.,  (4.13),  when 
there  are  no  simultaneous  observations,  and  block  lower  triangular  in  general,  e.g., 
(4.15).  Consider  the  specific  case  shown  in  Fig.  4.2.  By  writing  (4.41)  as  E'  =  EAq  1 , 
the  matrices  have  the  form 
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(4.42) 

It  can  be  seen  that  (4.42)  is  similar  in  form  to  M  =  UL,  where  U  and  L  are  upper 
and  lower  triangular  matrices,  respectively.  Thus  if  E'  can  be  factored  into  E  and 
Ag  l,  then  the  solution  to  the  multirate  Normal  equations  has  been  found. 

Most  linear  algebra  books  (e.g.,  [Ref.  65])  contain  LU  factorization  algorithms 
for  solving  the  decomposition  of  a  matrix  into  a  lower-upper  triangular  scheme  (i.  e., 


82 


M  =  LU).  For  the  purpose  of  this  dissertation,  however,  it  is  desired  to  factor  E' 
into  an  upper-lower  scheme.  In  addition,  this  factorization  must  be  general  enough  to 
allow  for  block  elimination  in  the  algorithm  (this  accounts  for  joint  observations  in  the 
multirate  system).  Thus  a  generalized  UL  factorization  algorithm  has  been  written 
and  provided  in  Appendix  E.  Using  this  generalized  UL  factorization  algorithm  to 
find  A0,  the  multirate  Normal  equations  can  be  solved.  A  significant  advantage  of 
decomposing  E'  into  E  and  Ag  1  is  that  the  only  matrix  that  requires  inversion  is  Ag  1. 
Since,  Ag  1  is  lower  triangular,  Ag  can  be  found  efficiently  using  forward  elimination 
which  does  not  explicitly  perform  a  full  matrix  inversion  [Ref.  66] . 

The  steps  needed  to  solve  for  the  multirate  Normal  equations  are  summarized 
in  Table  4.1. 


Table  4.1.  Steps  in  Solving  the  Multirate  Normal  Equations 

1.  Solve  (4.39)  using  the  multichannel  Levinson  recursion  to  find  E',  Ax, . . . ,  A q 

2.  Factor  E'  using  the  generalized  LIL  factorization  algorithm 
to  obtain  E  and  Ag  1 

3.  Invert  Ag  1  to  obtain  Ag 

4.  Compute  A,  =  A]  A0  for  i  —  1, . . . ,  Q 
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D.  SUMMARY 


In  this  chapter,  the  multirate  linear  prediction  equations  were  developed.  Also, 
an  efficient  method  for  calculating  the  linear  prediction  coefficients  using  the  multi¬ 
channel  Levinson  recursion  and  LU  factorization  was  propsed.  These  techniques  for 
calculating  the  linear  prediction  coefficients  are  used  in  the  next  chapter  on  multirate 
sequential  classification. 
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V.  CLASSIFICATION 


The  problem  of  distinguishing  between  two  or  more  different  types  of  signals 
(arising  from  different  sources)  is  known  as  classification.  The  classification  is  said  to 
be  “sequential”  when  the  signal  to  be  classified  is  observed  one  discrete  sample  at  a 
time.  After  a  new  sample  is  observed,  an  attempt  is  made  to  classify  the  signal  using 
all  of  the  samples  available  up  to  that  time.  If  a  classification  cannot  be  made,  then 
another  sample  is  taken  and  the  procedure  continues.  The  goal  is  to  classify  the  signal 
with  a  prescribed  probability  of  error  using  as  few  samples  as  possible.  The  general 
theory  of  sequential  classification  was  originally  and  extensively  developed  by  Wald 
in  1947  [Ref.  20].  It  has  since  been  applied  specifically  to  statistical  signal  processing 
and  described  for  the  single-channel  and  multichannel  cases  (e.g.,  [Ref.  22]  and  [Ref. 
21].)  The  focus  of  this  chapter  is  to  demonstrate  the  feasibility  of  developing  an 
algorithm  that  allows  for  multirate  observations  to  be  used  in  a  sequential  classifier. 

A.  SEQUENTIAL  CLASSIFICATION 

Before  developing  the  sequential  classifier  for  the  multirate  system,  it  is  nec¬ 
essary  to  provide  a  brief  description  of  the  sequential  classification  process.  Conven¬ 
tional  methods  of  classification  usually  involve  a  fixed  block  of  data.  This  method 
is  known  as  classifying  with  a  fixed  sample  size.  Statistically  optimal  methods  then 
involve  developing  the  ratio  of  likelihood  functions  for  the  two  classes  and  comparing 
that  likelihood  ratio  to  a  fixed  threshold.  This  threshold  is  chosen  to  optimize  some 
criterion,  such  as  total  probability  of  error  (Bayes  test)  or  maximizing  the  probability 
of  correctly  choosing  one  class  while  holding  the  probability  of  error  on  the  other 
class  fixed  (Neyman- Pearson  test)  [Ref.  22,  67,  68].  Wald  developed  a  method,  called 
the  sequential  probability  ratio  test  (SPRT),  which  allowed  one  to  successively  add 
samples  to  data  set  being  used  to  perform  the  classification.  If  classification  is  not 
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possible  at  that  time,  another  sample  is  added.  This  process  is  repeated  until  a  clas¬ 
sification  can  be  made  with  some  fixed  desired  probability  of  error.  For  an  extensive 
discussion  on  sequential  analysis,  see  [Ref.  20]. 

To  begin  the  discussion  of  sequential  classification,  assume  that  the  signal  of 
interest,  x ,  belongs  to  one  of  only  two  classes,  class  0  or  class  1.  Let  H0  be  the 
hypothesis  that  the  signal  belongs  to  class  0  and  Hi  the  hypothesis  that  the  signal 
belongs  to  class  1.  Let  the  vector  of  observations  be  denoted  by 


x„  = 


Xr 


(5.1) 


where  Xj  is  the  set  of  observations  of  x  available  at  time  j,  and  let  /(xn|iL,)  be  the 
probability  density  function  for  xn  given  that  the  observations  are  from  class  i.  In 
the  Wald  SPRT,  the  likelihood  ratio  is  defined  as 


^  /(xw|iLi) 
/(Xn|i?o) 


(5.2) 


and  two  positive  constants  or  thresholds  A  and  B  (with  A  >  B)  are  chosen  to  classify 
the  signal  xn  as  either  class  0  or  class  1  with  some  desired  probabilities  of  error.  The 


test  performed  is  as  follows:  if 

/(Xn|ffl)  >  A 

f(xn\H0)  ~ 

then  xn  is  classified  as  class  1;  while  if 


(5.3) 


/(*■■  l-ffi)  „ 
/(x„|iJo)  ~ 


(5.4) 


then  xn  is  classified  as  class  0.  If  the  ratio  of  the  conditional  probabilities  falls  between 


A  and  B ,  that  is,  if 


B  < 


f  (xw  |  H\ ) 
f  (Xn  |  Hq  ) 


<  A, 


(5.5) 


then  another  sample  is  taken  and  the  process  is  repeated  with 


xn+l  = 

—n 

— n+1 

Proper  selection  of  the  values  for  A  and  B  is  discussed  below. 

In  many  cases,  it  is  more  convenient  to  work  with  the  logarithm  of  the  like¬ 
lihood  ratio  rather  than  the  likelihood  ratio  itself.  Since  the  logarithm  is  a  strictly 
monotonic  function,  this  does  not  change  the  essential  nature  of  the  SPRT.  In  the 
particular  case  that  follows,  it  is  most  convenient  to  deal  with  the  negative  logarithm 
and  define 

/in(Xn)  =  -2  ln(/n(x„))  (5.6) 

with 


ta  =  -2  ln(,4) 

(5.7) 

tb  =  -2  ln(R) 

(5.8) 

The  inequality  then  reverses  and  the  SPRT  decision  becomes 

<  ta  choose  Hi 

Mxn) 

>  tb  choose  H0 


where 


^n(^n)  —  2  111 


fxn\Hi  (xn|-^l) 
/xn  |  Ho  ( xn  |  Hq  ) 


Otherwise  take  another  observation. 


(5.9) 

(5.10) 


The  choices  for  the  thresholds  A  and  B  are  related  to  the  probabilities  of  error 
for  class  0  and  class  1.  Using  the  method  described  in  [Ref.  22,  21],  the  relations 
between  thresholds  and  probabilities  of  error  can  be  derived  as  follows.  Let  eo  be  the 
probability  of  a  decision  error  given  the  data  is  from  class  0  and  e\  be  the  probability 
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of  error  given  the  data  is  from  class  1.  Assume  that  the  upper  threshold  is  crossed 
after  n  observations;  that  is,  at  time  n  the  ratio  /„(xn)  defined  by  (5.2)  is  greater 
than  but  approximately  equal  to  A: 


^n(Xn)  >  A 

Substituting  (5.2)  into  (5.11)  and  rearranging  yields 

f  (xn | H\ )  >  Af(xn\H0). 


(5.11) 


(5.12) 


The  condition  (5.12)  defines  a  region  of  the  xn  space  (Ri)  where  the  data  is  classified 
as  class  1  (Hi).  Integrating  both  sides  over  this  region  yields 


/(xn|/A)dxn  >  A 


f(xn\H0)dxn 


(5.13) 


Now,  since  R\  is  taken  to  be  the  region  of  measurement  space  where  xn  is  classified 
as  class  1,  the  left  hand  side  of  (5.13)  is  the  probability  of  correctly  classifying  class 
1,  namely 


f('x-n\Hi)d'X.n 


1  —  €\. 


(5.14) 


The  integral  on  the  right  hand  side  of  (5.13)  is  then  the  probability  of  incorrectly 
classifying  class  0, 


f(xn\H0)dxn 


eo 


(5.15) 


Substituting  (5.14)  and  (5.15)  into  (5.13)  yields 


1  —  ei  ^  Aeo 


(5.16) 


or 

A  =  (5.17) 

eo 

where  the  inequality  has  been  replaced  by  an  equal  sign  since  the  threshold  A  is  most 
restrictive  at  that  value.  Similarly,  assuming 


l„(x„)  <  B 


(5.18) 


and  integrating  over  the  region  where  xn  is  classihed  as  class  0  yields 

B  =  (5.19) 

I  —  e0 

Equations  (5.17)  and  (5.19)  provide  the  necessary  relations  to  allow  one  to  choose 
appropriate  values  of  thresholds  A  and  B  based  on  desired  probabilities  of  error. 

B.  ALGORITHM  DEVELOPMENT 

When  the  observations  in  the  SPRT  are  jointly  Gaussian  for  both  classes, 
then  a  convenient  recursive  classification  algorithm  can  be  developed  that  involves 
linear  prediction.  This  recursive  form  has  been  developed  in  [Ref.  22,  21]  for  the 
case  of  a  single  observation  sequence,  which  is  scalar-valued  or  vector-valued  (the 
multichannel  case).  Here,  the  method  is  extended  to  the  case  of  multiple  channels 
sampled  at  different  rates.  Thus,  at  any  given  epoch  on  the  system  time  scale,  there 
may  be  one  or  more  simultaneous  observations  available  from  the  various  channels 
or  possibly  no  observations.  The  algorithm  to  be  described  takes  this  more  general 
situation  into  account. 

The  development  of  the  recursive  form  of  the  quadratic  classifier  for  the  mul¬ 
tirate  system  closely  follows  the  development  for  the  single-channel  and  multichannel 
cases  developed  in  [Ref.  22]  and  [Ref.  21].  It  is  necessary,  however,  to  make  modifi¬ 
cations  to  some  of  the  parameters,  to  account  for  the  periodic  nature  of  the  multirate 
system.  Thus,  some  parameters  that  are  constant  in  the  single-channel  and  multi¬ 
channel  cases  become  periodically  time- varying  variables  in  the  multirate  case. 

To  start  the  development,  assume  that  the  signals  associated  with  both  classes 
are  jointly  Gaussian  and  that  the  observations  for  all  channels  are  collected  and 
ordered  such  that  xn  contains  all  observations  from  time  0  to  n.  In  addition,  the 
observation  vector  (5.1)  can  be  written  recursively  as 

(5.20) 


x„  = 


1 

(k) 
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where  the  superscript  k  is  used  to  indicate  that  the  size  of  the  observation  vector 
r?  is  a  periodically  time-varying  parameter  since  the  number  of  observed  channels 
at  any  time  step  varies  periodically.  If  the  size  of  xi!^  (the  number  of  observations  at 
time  n )  is  defined  as  Sn ^  and  the  total  number  of  observations  from  time  0  to  n  —  1 
is  defined  as  Nn- 1,  then  the  total  number  of  observations  at  time  n  can  be  expressed 
as 

Nn  =  Nr^1+Sjlkl  (5.21) 

As  an  example,  consider  the  two-channel  system  where  the  first  channel  is 
decimated  by  a  factor  of  2  and  the  second  channel  is  decimated  by  a  factor  of  3,  as 
shown  in  Fig.  5.1.  The  system  period  as  defined  in  Chapter  11  is  K  —  6.  Thus,  since 
n=k  (mod  Ii),  the  system  phase,  k,  varies  periodically  between  0  and  5. 
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Figure  5.1.  Two-channel  Multirate  System. 


This  can  be  illustrated  more  specifically  as  follows  (refer  to  Fig.  5.1).  Assum¬ 
ing  that  observations  through  n  —  29  have  been  collected,  the  associated  observation 
vectors  for  time  steps  n  =  30,  31, . . . ,  35  would  be  denoted  by 


X29 

X30 

X31 

X30  = 

Mo) 

X31  = 

„(!) 

X32  = 

_(2) 

T30 

T31 

i32 

X32 

X33 

X34 

X33  = 

M3) 

X34  = 

_(4) 

X35  = 

„(5) 

—33 

—34 

i35 
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where  the  vectors  of  observations  at  these  times  are 


T(o) 

—30  — 


aii  [30] 
x2[30] 


T,3i)  —  empty  X32  ~  x\  [32] 


^33  =^[33]  =  empty 

The  corresponding  number  of  observations,  S ^ ,  at  these  times  are 


c(0)  _  C\  rr(l)  _  n  c(2)  _  1 

°30  —  z  °31  —  u  °32  —  L 


c(3)  —  1  c(4)  —  1  o(5)  _  n 

°33  _  1  °34  —  1  °35  —  u 


For  any  given  time  ln\  the  probability  density  function  for  the  collection  of 
observations  is  given  by  the  multirate  Gaussian  form 


/x„(xn)  — 


(27r)JV"/2|Cn|1/2 


p-i(x„-m„)TC„1(x„-m„) 


(5.22) 


where  mn  =  £{xn}  is  the  mean  vector  and  Cn  =  £{(xn  —  mn)(xn  —  mn)T}  is  the 
covariance  matrix  for  the  observations.  The  mean  vector  and  covariance  matrix  can 


then  be  written  as 


mn_i 


(5.23) 


C„  = 


p  T> 

'^'n—  1 

-p(Al)T  yi(fc) 


(5.24) 


where  the  partitioning  corresponds  to  the  partitioning  defined  in  (5.20).  By  using  a 
formula  for  matrix  partitioning  [Ref.  69],  the  inverse  covariance  matrix,  C"1,  can  be 
represented  in  terms  of  the  inverse  covariance  matrix,  [ ,  as 


c.-1  = 


I  — G.Lfc)  c~\  0 


I  0 


0T  I 


r»T  1  /~((fc)T  T 

U  £jn  vjn  1 


(5.25) 
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or  equivalently  as 


c-1,  0 

+ 

(k) 

jy(fc)  1 
n 

-G  ^  I 

(5.26) 

0T  0 

I 

L  J 

where  G.^  and  are  defined  as 


C 


-i 
n—  1 


y(k)  _  f{(fc)Tn-1  ■R(fc) 

■L~‘n  ^n—  lAl/n  ■ 


(5.27) 

(5.28) 


Since  the  determinants  of  the  upper  and  lower  triangular  matrices  of  (5.25)  are  unity, 
the  determinant  of  C,”1  is  given  by 


ic-1!  =  ic 


-1 
n—  1 


|Eifc)' 


or 


(5.29) 


The  partitionings  of  (5.1)  and  (5.23)  through  (5.29)  provide  all  the  necessary  relations 
for  the  recursive  computation  of  the  density  function  (5.22)  at  each  time  step. 

Now  that  the  probability  density  function  for  the  multirate  system  has  been  de¬ 
scribed,  the  probability  density  function  for  a  given  class  within  the  multirate  system 
can  be  discussed.  Consider  observations  from  a  known  class  i  with  a  corresponding 
mean  vector 


rn 


(0 


rn 


m 


(0 

n—1 

«(*) 


(5.30) 


and  covariance  matrix 


^n- 1 

■p  w  w 

T3  (*)(^)T 

-K-n 

yWW 

(5.31) 


where  the  superscript  (i)  denotes  association  with  class  i  and  the  partitioning  corre¬ 
sponds  to  the  partitioning  defined  in  (5.23)  and  (5.24).  The  vector  of  observations 
with  the  mean  removed  can  be  denoted  as 


xh)  —  x 

-*~n  ■A-n 


—  m,il) 


(5.32) 
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or 


xW  = 


xw 

xn—  1 

_(0(fc) 

(i) 

Xji-i  - 

r(U  _  ra(')(fc) 

djji  LLLn 


(5.33) 


By  using  the  form  (5.22)  and  substituting  the  relations  (5.21),  (5.32),  (5.33),  (5.26) 
and  (5.29),  the  probability  density  function  for  class  i  at  time  n  can  be  written  as 


cl'1  (xn-mj,'1] 


(27T)(^)/2|Cii)|1/2 


1  fx(i)Tr(i)_1x(i)  'l 

3  2  \Xn— 1  y~'n—l  Xn-1  ) 


X 


(27r)ivn-1/2|C«i|i/2  (27r)s™fc)/2|E«(fc)|i/2 

3_1  (2«w_gwtx«  ^  (534) 


x  e 


The  first  term  on  the  right  represents  /Xn_1|f/-i  while  the  second  term  on  the  right 
represents  a  conditional  probability  density  function  ,/x  |xn_i/7.r  In  other  words,  (5.34) 
can  be  written  as 

/x„|Hi  =  fxn-i\Hi  '  f  =  1, 2  (5.35) 

where 

„  1  _i  fJOTptrhw  1 

/xn_i|iJi  —  77  777 - „,„e  2U-1  »-i  n-i)  (5.36) 

and 

1 


J2n\ 


xn—  lHi 


(27r)sifc)/2|E®(fc)|1/2 


(27r)^n-l/2|C«1|1/2 

lf_(0(fc)  rWT  <%)  lTp(i)(fc)_1  /'_(*)(&)  r(k)T  (i)  \ 

g  2  71  Xro  —  1  ^  “n  t— n  Xn-1  / 


Using  these  results,  the  statistic  /rn(xn)  defined  in  (5.10)  becomes 
Mx„)  =  -2111  [  X—  ift(x-i|gi) . 

/x,1_i|H0(xn-l|h/o)  Jxn|xn_ii?0 


_  2  In  ^ Xn~  1  l^i  ^ x 77  1 1  ^^1 )  2  in  ^ ~n  lXn~ 1  ^1 

/xn_i|ff0(Xn-l|^o)  /xn|xn_ii?0 


(5.37) 


(5.38) 
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This  can  be  written  as 


where 


1  +  A/in(x^}|xn_i), 

(5.39) 

Q  i  ^  —n  lxn— 1  -^1 

/iJXn-lHo 

(5.40) 

Substituting  (5.37)  into  (5.40)  results  in  the  expression 


1  4  fe"*  - 

1  -vr  .  1  - 

r(fc)(l): 

rx(l)  42  /  (0) 

xn— 1 / 

n M(0)‘ 

'Jn 

rx(°)  j2 

Xn~l)  ,  i 

E^')(1)| 

|  Xn—  1 )  — 

ipW(i) 

Enfc)(0) 

HI 

Enfc)(0)] 

(5.41) 


Since  the  first  observation  vector  is  labeled  x0  and  the  statistic  ho(x0)  = 
h- i(x_i)  +  A/i0(a;o0)|x-i)>  the  SPRT  is  first  initialized  with  X-\  =  0  and  the  statistic 
h- i(x_i)  is  empty.  Thereafter,  hn  is  computed  from  (5.39)  and  (5.41)  recursively. 
This  statistic  is  used  in  the  decision  rule  (5.9)  to  perform  the  classification. 


C.  CLASSIFICATION  USING  LINEAR  PREDICTION 


The  equations  for  the  multirate  sequential  classifier  were  developed  assuming 
that  the  observations  were  jointly  Gaussian.  The  sequential  classifier,  however,  can 
also  be  interpreted  in  terms  of  linear  prediction.  By  rearranging  (5.27)  and  (5.28) 
into 

R‘‘)-C,.1gW  =0  (5.42) 

=  (5.43) 


these  two  equations  can  be  combined  in  matrix  form  to  write 


1 

h 

i 

1 

H 

^2 

'pH 
_ 1 

1 

1 

e 

u 

i 

2e 

0 

1 

- 1 

o 

(5.44) 


These  equations  bear  a  striking  similarity  to  the  Normal  equations  of  linear 
prediction  (see  (4.2)).  Indeed,  the  matrix  -G®  can  be  interpreted  as  the  prediction 
filter  coefficients  for  predicting  the  new  observations,  and  is  the  prediction  error 
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covariance  matrix.  More  explicitly  and  relate  to  the  parameters  in  (4.9) 


by 


and 


r  “i 

I  1 

0  I  . 

■  1  0 

Ao 

Ax 

I 

-gLa_1)  1 

1 

I  1  . 

.  I  0 

= 

- G (?"2>  I  . 

■  1  0 

.  1  I 

.aq. 

I 

1  r(0) 

E^_1) 

-L'n 

x 

.  .  1  X 

0 

l  i 

C-in  \ 

. .  |  X 

■ .  I  x 

0 

1  0  1  . 

1  E® 

(5.45) 


(5.46) 


The  block  form  in  (5.45)  and  (5.46)  may  be  a  little  misleading.  In  most  cases, 
the  width  of  the  blocks  is  just  one  column  and  the  becomes  a  “1.”  The  matrices 
in  (5.45)  and  (5.46)  have  been  represented  in  block  form  to  allow  for  possible  multiple 
observations  at  each  time  step.  The  parameters  -Gn'  and  thus  represent  the 
periodically  time-varying  parameters  used  to  predict  each  new  observation  (or  set  of 
observations)  within  a  block. 

When  a  sufficiently  high  order  Q  is  chosen  for  the  prediction,  the  parameters 
Gn'1  and  E,^  do  not  depend  explicitly  on  the  time  variable  “n”  but  only  on  the  phase 
as  indicated  by  the  superscript  k. 


D.  ALGORITHM  SUMMARY 

This  section  summarizes  the  multirate  sequential  classifier  algorithm  for  ease 
of  implementation.  The  recursive  classification  algorithm  can  be  divided  into  two 
parts:  1)  training  and  2)  testing. 

In  the  first  part,  the  system  is  trained  to  distinguish  between  two  classes  using 
a  priori  information.  Using  known  class  training  data,  the  classifier  parameters  (mean 
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vectors,  prediction  coefficient  matrices  and  prediction  error  variance  matrices)  for  each 
class  are  estimated.  In  addition,  the  thresholds  used  in  classification  are  determined 
and  adjusted  experimentally.  These  parameters  are  summarized  in  Table  5.1. 


Table  5.1.  Summary  of  Classifier  Training  Parameters 


Class  i  Parameters: 

mean  vectors,  mP1  =  £i[xn*^] 

i  =  0,1 

Prediction  coefficients,  Gn^fc) 

Prediction  Error  Variances, 

Classification  Thresholds: 

A  (Class  1)  and  B  (Class  0) 

In  the  second  part  new  (test)  data  is  presented  to  the  classifier.  The  associated 
class  errors  (en^  and  )  for  the  given  signal,  xn,  are  calculated  and  used 

along  with  the  class  prediction  error  variances  to  calculate  the  SPRT  statistic,  hn , 
for  comparison  with  the  class  thresholds.  When  the  SPRT  statistic  crosses  a  class 
threshold,  the  classifier  declares  the  classification;  otherwise,  the  process  is  repeated. 
These  steps  are  shown  in  Table  5.2. 


Table  5.2.  Steps  in  Sequential  Classification 
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The  performance  of  the  classifier  can  be  evaluated  with  two  key  measures  of 
performance:  the  number  of  correct  classifications  and  the  length  of  data  (length  of 
time)  required  to  make  a  classification.  By  varying  different  properties  of  the  classifier 
and  environment  (e.g.,  signal-to-noise  ratio),  the  effects  of  these  different  properties  on 
the  classifier  can  be  evaluated  using  the  two  performance  measures.  The  remainder 
of  this  chapter  describes  the  simulations  and  results  for  varying  conditions  of  the 
classifier  and  environment. 

E.  SIMULATION  SETUP  AND  PARAMETERS 

In  order  to  test  the  proof  of  concept  for  the  multirate  sequential  classifier 
multiple  simulations  were  conducted.  The  test  data  consisted  of  recordings  of  sound 
from  one  propellor  plane  and  three  different  A-10  jet  aircraft,  with  noise  added.  The 
parameters  that  were  varied  during  the  simulations  were  the  signal-to-noise  ratio 
(SNR)  associated  with  the  training  and  target  data  sequences,  and  the  length  of 
the  training  sequences  used  to  develop  the  classifier  parameters.  Simulations  were 
conducted  to  compare  the  single-channel,  multichannel  and  multirate  cases. 

1.  Test  Data  Characteristics 

This  subsection  describes  the  characteristics  of  the  test  data.  The  original 
audio  data  sequences  were  recorded  at  44,100  Hz.  The  data  was  initially  low-pass 
filtered  and  downsampled  by  a  factor  of  10  to  achieve  a  sampling  rate  of  4,410  Hz. 
(The  observed  data  for  the  simulation  is  then  further  decimated.)  A  spectral  estimate 
was  then  computed  for  each  of  the  four  data  sequences  prior  to  adding  any  additional 
noise,  in  order  to  determine  any  unique  characteristics.  Figure  5.2  shows  the  spectral 
plots  of  each  sequence  and  Table  5.3  summarizes  the  dominant  frequencies,  and  a 
secondary  frequency  if  noticeable,  that  were  measured  from  the  spectral  plots. 

From  Fig.  5.2,  it  can  be  seen  that  the  propellor  plane  has  secondary  frequency 
at  a  lower  frequency  than  the  peak  frequency  of  192.0  Hz.  A  secondary  frequency  is 
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Frequency  (Hz) 
(d)  A-IOJet  No.  3 


Figure  5.2.  Spectral  Plots 


Table  5.3.  Spectral  Information  for  Data  Sequences 


Propellor  Plane 

A-10  Jet  No.  1 

A-10  Jet  No.  2 

A-10  Jet  No.  3 

Dominant 
Freq  (Hz) 

192.0 

145.0 

82.5 

84.6 

Secondary 
Freq  (Hz) 

106.9 

N/A 

N/A 

210.4 
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also  seen  on  the  A-10  Jet  No.  3  at  210.4  Hz.  The  A-10  Jet  No.  1  has  many  peaks  close 
to  each  other.  However,  the  magnitude  of  these  peaks  diminish  as  the  distance  from 
the  main  peak  increases.  Thus  they  have  not  been  considered  as  secondary  peaks. 

In  addition  to  the  spectral  estimate  of  the  test  data,  the  autocorrelation  of  each 
data  sequence  was  calculated  to  examine  any  significant  differences  in  correlation. 
The  autocorrelation  plots  were  used  to  assist  in  selecting  appropriate  filter  orders  for 
calculating  the  prediction  orders  used  in  the  classifier.  One  can  see  from  Fig.  5.3  that 
the  three  A-10  Jet  planes  have  similar,  but  not  identical,  correlation  structure  within 
the  first  ten  samples  (dashed  lines  in  Fig.  5.3.)  This  similar  structure  led  to  choosing 
a  prediction  filter  order  larger  than  ten,  as  discussed  in  the  next  subsection. 


Lag  Lag 

(a)  Autocorrelation  for  Propeller  Plane  (b)  Autocorrelation  for  A-1 0  Jet  No.  1 


Lag  Lag 

(c)  Autocorrelation  for  A-10  Jet  No.  2  (d)  Autocorrelation  for  A-1 0  Jet  No.  3 


Figure  5.3.  Estimated  Autocorrelation  Function  for  Aircraft  Data  at  4,410  Hz  sam¬ 
pling  rate.  Dashed  lines  show  ±10  lag  values. 
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2.  Simulation  Models 

For  the  simulations  conducted,  three  scenarios  were  considered.  These  scenar¬ 
ios  considered  classification  using  a  single-channel  system,  the  multichannel  system 
and  the  multirate  system.  All  scenarios  began  with  the  recorded  data  sampled  at 
4410  Hz.  For  the  single-channel  scenario,  the  recorded  data  was  decimated  by  a  fac¬ 
tor  of  two,  resulting  in  a  sampling  rate  of  2205  Hz.  This  was  done  to  ensure  that 
the  highest  rate  in  all  three  scenarios  would  be  the  same.  Additive  white  Gaussian 
noise  with  a  mean  of  zero  and  a  variance  of  a ^  was  then  added  to  achieve  a  desired 
Signal-to-Noise  ratio  set  by  the  simulation  parameters.  The  proper  selection  of  the 
value  of  cbj  is  described  below.  The  single-channel  scenario  is  shown  in  Fig.  5.4. 


Channel  1 
(2205  Hz) 


Figure  5.4.  Observation  Model  for  the  Single  Channel  Scenario. 


The  simulation  model  for  the  multichannel  scenario  is  similar  to  the  single¬ 
channel  scenario,  with  the  exception  that  after  the  recorded  data  is  decimated  by  a 
factor  of  two,  the  signal  is  split  into  two  channels.  Noise  was  then  added  to  each 
channel  to  achieve  a  desired  SNR.  The  model  for  the  multichannel  scenario  is  shown 
in  Fig.  5.5. 

For  the  multirate  scenario,  Fig.  5.6,  the  high  rate  channel  was  decimated  by  a 
factor  of  2  and  the  low  rate  channel  was  decimated  by  a  factor  of  3.  Noise  was  then 
added  to  each  channel  to  achieve  a  desired  signal-to-noise  ratio. 

In  order  to  produce  the  desired  SNR  for  these  models,  the  noise  power  must 
be  calculated  using  the  desired  SNR  in  decibels  and  the  estimated  power  of  the 
appropriate  channel.  If  the  output  of  the  decimator  is  s[mx],  where  mx  is  the  time 
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Noise  1 


Channel  1 
(2205  Hz) 


Channel  2 
(2205  Hz) 


Noise  2 


Figure  5.5.  Observation  Model  for  the  Multichannel  Scenario. 


Noise  1 


Channel  1 
(2205  Hz) 


Channel  2 
(1470  Hz) 


Noise  2 


Figure  5.6.  Observation  Model  for  the  Multirate  Scenario. 
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index  for  the  decimated  signal,  r)[mx]  is  the  noise  signal,  and  the  observed  signal  x[mx] 

is 

x[mx]  =  s[mx\  +  r)[mx ], 

then  the  noise  power,  a can  be  calculated  as 

mx  —  i 

£  E 

a 2  = _ _ 

V  1Q-SNR/10 

Once  the  noise  power  has  been  calculated,  this  value  can  be  used  to  generate  the 
appropriate  noise  signal,  rj[mx]. 

3.  Prediction  Coefficients 

When  calculating  the  prediction  coefficients  for  the  training  data,  an  appro¬ 
priate  filter  order  has  to  be  chosen.  In  general,  the  filter  order  should  be  large  enough 
to  exploit  differences  in  correlation  between  the  two  classes  of  signals  being  tested. 
It  should  be  small  enough,  however,  so  that  it  does  not  add  excessive  computational 
burden  to  the  simulation  and  does  not  become  sensitive  to  small  artifacts  of  the  train¬ 
ing  data.  In  addition,  as  the  filter  order  increases,  the  determinant  of  the  prediction 
error  covariance,  which  measures  the  quality  of  the  prediction,  typically  approaches 
a  value  such  that  any  further  increases  in  filter  order  do  not  result  in  an  appreciable 
improvement  in  performance. 

From  Fig.  5.3,  it  can  be  seen  that  the  propellor  plane  is  uniquely  different  in  its 
correlation  structure,  but  the  three  A-10  jets  have  similar  structures  for  lags  smaller 
then  ten.  Thus  choosing  an  order  for  the  linear  prediction  filter  less  than  ten  would 
probably  not  provide  enough  uniqueness  in  the  three  A-10  jet  parameters  to  be  able 
to  adequately  classify  the  different  planes.  However,  the  correlation  shapes  at  lags 
larger  than  ten  are  sufficiently  dissimilar  to  exploit  these  differences  in  correlations. 
Although  larger  filter  orders  could  have  been  used,  the  filter  order  was  chosen  to  be 
40  to  minimize  excess  computations. 
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The  prediction  coefficients  and  prediction  error  variances  were  calculated  in  the 
manner  described  in  Chapter  IV,  using  the  multichannel  Levinson  recursion.  Since 
the  decimation  rates  for  the  two  channels  in  the  multirate  scenario  are  two  and  three, 
the  system  period,  K ,  is  equal  to  6.  In  addition,  there  are  a  total  of  five  observations 
that  occur  within  one  system  period,  three  for  the  high-rate  channel  and  two  for  the 
low-rate  channel.  The  sampling  pattern  is  shown  in  Fig.  5.7.  It  is  assumed  the  data 
is  aligned  so  that  at  the  initial  observation,  both  channels  produce  observations. 


System  Period  K  =  6 

r  ^  \ 


0  12345678 


Channel  1 


Channel  2 


O— 0—0— 0—0 

Kx  =2 


K2  =  3  time 


Figure  5.7.  Illustration  of  Decimation  for  the  Multirate  Scenario. 


4.  Training  Data 

When  calculating  the  prediction  coefficients,  a  set  of  training  data  for  each 
class  had  to  be  used.  This  data  consisted  of  a  portion  of  the  total  data  signal  used 
for  each  class,  as  shown  shaded  in  Fig.  5.8.  The  length  of  this  portion  was  varied 
according  to  the  scenario  parameters.  The  rest  of  the  signal  was  then  used  for  the 
testing  portion  of  the  simulation. 

During  the  testing  phase,  the  simulation  starts  at  the  beginning  of  the  test 
portion  of  the  input  signal  and  continues  to  take  observations  until  a  classification  is 
made.  This  is  considered  to  be  one  “trial.”  A  second  trial  is  then  begun  starting  with 
the  observation  corresponding  to  the  data  point  just  after  the  data  point  at  which 
the  first  classification  was  made,  as  shown  in  Fig.  5.8.  This  process  is  continued 
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until  50  trials  are  conducted.  The  length  of  the  data  was  sufficient  to  ensure  that  the 
simulation  was  able  to  conduct  50  trials  before  reaching  the  end  of  the  data. 


Start  of 
Classifier  Data 

Figure  5.8.  Training  Data. 


Although  the  sequential  classifier  can  in  theory  classify  data  based  on  as  little 
as  one  single  observation,  in  practice  it  is  better  to  allow  some  minimum  number  of 
observations  to  be  collected  before  a  classification  is  attempted.  In  these  experiments, 
this  minimum  data  length  was  linked  to  the  length  of  the  linear  prediction  filter.  In 
particular,  the  algorithm  was  designed  not  to  classify  the  target  until  the  number  of 
observations  exceeded  the  filter  order  by  one  observation.  It  was  decided  to  implement 
the  algorithm  this  way  to  avoid  any  errant  classifications  that  might  occur  while  the 
algorithm  collects  the  observations  needed  to  make  the  predictions  at  full  order. 


F.  SIMULATION  RESULTS  AND  ANALYSIS 

The  results  for  the  multirate  scenarios  are  presented  in  the  following  subsec¬ 
tion;  the  key  parameters  are  analyzed  in  the  last  subsection. 
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1. 


Multirate  Scenario  Results 


Since  there  were  four  classes  of  data  (one  propcllor  and  three  jet  planes)  the 
classifier  was  tested  using  two  classes  at  a  time.  This  resulted  in  six  sets  of  experi¬ 
ments,  whose  results  are  listed  in  Tables  5.4  through  5.9  of  this  section. 

In  order  to  make  comparisons,  specific  parameters  were  varied  for  the  simula¬ 
tions.  These  were  the  signal-to-noise  ratio  (SNR)  of  the  target  data  and  the  length 
of  data  sequences  used  to  calculate  class  parameters.  In  all  cases,  the  training  data 
was  assumed  to  be  noiseless  but  white  noise  was  added  to  the  test  data  to  achieve  a 
desired  SNR. 

Table  5.4.  Classification  Results:  Propellor  vs.  A-10  Jet  No.  1 


SNR 

Data 

Training 

Length 

Propellor 

No.  of  Correct 

Classifications 
(out  of  50) 

Jet  No.  1 

No.  of  Correct 

Classifications 
(out  of  50) 

Average  Time 
to  Classify 
Propellor 
(msec) 

Average  Time 
to  Classify 

Jet  No.  1 
(msec) 

40  dB 

6250 

50 

50 

9.52 

9.52 

2500 

50 

50 

9.52 

9.52 

20  dB 

6250 

50 

50 

9.52 

9.59 

2500 

50 

39 

9.52 

9.52 

10  dB 

6250 

50 

21 

9.52 

9.52 

2500 

50 

10 

9.52 

9.52 

Table  5.4  shows  the  classification  results  for  the  propellor  plane  versus  the 
A-10  Jet  No.  1.  This  table  lists  both  the  number  of  correct  classifications  for  each 
target  as  well  as  the  average  time  to  classify  the  target.  The  average  time  of  9.52  msec 
listed  in  the  table  corresponds  to  the  minimum  observation  condition  of  41  samples 
(see  Subsection  V.E.4).  It  can  be  seen  from  Table  5.4  that  the  classifier  makes  its 
decision  almost  immediately  after  reaching  the  full  filter  order  under  all  cases.  Thus, 
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for  the  case  of  the  propellor  plane  versus  Jet  No.  1,  the  advantage  of  a  sequential 
classifier  is  not  realized.  This  is  not  true  for  all  the  test  scenarios. 

From  Table  5.4,  one  can  see  that  the  algorithm  is  able  to  classify  the  propellor 
under  all  conditions,  and  almost  immediately  after  collecting  enough  data  to  reach 
full  order.  In  an  environment  with  an  SNR  of  40  dB,  the  classifier  is  able  to  classify 
the  propellor  plane  nearly  perfectly.  As  the  SNR  is  reduced  to  20  dB,  the  success 
rate  reduces  to  a  little  below  80%  for  the  short  training  data  length  and  under  15% 
for  the  high  training  data  length.  With  the  middle  length  training  data  classification 
is  perfect.  Proportionately  similar  reductions  occur  as  the  SNR  is  dropped  further  to 
10  dB.  The  results  degrade  significantly  at  an  SNR  of  10  dB.  Namely,  the  classifier 
tends  to  classify  both  signals  as  a  propellor  plane;  this  puts  into  question  whether 
the  classifier  appropriately  classifies  the  propellor  plane  at  10  dB  or  defaults  to  the 
propellor  plane  because  of  the  signal  degradation.  In  almost  all  cases  classification 
occurred  with  the  minimum  allowable  number  of  observations.  For  SNR  below  5 
dB,  classification  degraded  severely.  Therefore,  the  results  are  not  presented  in  these 
tables. 


Table  5.5.  Classification  Results:  Propellor  vs.  A-10  Jet  No.  2 


SNR 

Data 

Training 

Length 

Propellor 

No.  of  Correct 

Classifications 
(out  of  50) 

A-10  No.  2 

No.  of  Correct 

Classifications 
(out  of  50) 

Average  Time 
to  Classify 
Propellor 
(msec) 

Average  Time 
to  Classify 
A-10  No.  2 
(msec) 

40  dB 

6250 

50 

50 

9.56 

9.52 

2500 

50 

50 

9.53 

9.52 

20  dB 

6250 

50 

36 

9.52 

9.53 

2500 

50 

8 

9.52 

9.53 

10  dB 

6250 

50 

2 

9.52 

9.52 

2500 

50 

0 

9.52 

9.52 

The  results  of  classification  for  the  propellor  plane  and  the  A-10  Jet  No.  2 
are  contained  in  Table  5.5.  Similar  to  the  testing  between  the  propellor  plane  and 
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the  A-10  Jet  No.  1,  the  algorithm  was  able  to  correctly  classify  the  propellor  plane 
nearly  perfectly  under  all  conditions  simulated.  In  addition,  the  classification  occurred 
almost  immediately  after  reaching  the  minimum  number  of  observations  allowed.  The 
number  of  correct  classifications  for  Jet  No.  2  was  slightly  less  than  that  for  Jet  No.  1. 
As  with  Jet  No.  1,  the  performance  of  Jet  No.  2  dropped  appreciably  as  the  SNR 
was  reduced  to  20  dB.  Under  the  conditions  with  an  SNR  of  10  dB,  the  classifier 
performed  poorly  when  trying  to  classify  Jet  No.  2. 


Table  5.6.  Classification  Results:  Propellor  vs.  A-10  Jet  No.  3 


SNR 

Data 

Training 

Length 

Propellor 

No.  of  Correct 

Classifications 
(out  of  50) 

A-10  No.  3 

No.  of  Correct 

Classifications 
(out  of  50) 

Average  Time 
to  Classify 
Propellor 
(msec) 

Average  Time 
to  Classify 
A-10  No.  3 
(msec) 

40  dB 

6250 

50 

50 

9.52 

9.52 

2500 

50 

50 

9.52 

9.52 

20  dB 

6250 

50 

45 

9.52 

9.53 

2500 

50 

35 

9.52 

9.53 

10  dB 

6250 

50 

2 

9.52 

9.52 

2500 

50 

0 

9.52 

9.52 

For  the  experiments  conducted  between  the  propellor  plane  and  the  A-10  Jet 
No.  3  shown  in  Table  5.6,  the  algorithm  once  again  was  able  to  classify  the  propellor 
plane  one  hundred  percent  of  the  time  and  with  the  minimum  number  of  observations 
allowed.  The  results  for  Jet  No.  3  were  somewhat  between  the  result  for  the  other 
two  A-10  jets.  For  the  simulations  with  an  SNR  of  40  dB,  the  algorithm  was  able 
to  classify  Jet  No.  3  perfectly.  As  the  SNR  was  reduced  to  20  dB,  the  performance 
dropped  to  between  70%  and  90%.  The  performance  then  significantly  dropped  to 
almost  no  correct  classifications  with  an  SNR  of  10  dB. 

Besides  comparing  the  jet  planes  to  the  propellor  plane,  it  was  also  of  interest 
to  see  if  the  classifier  could  discriminate  between  the  different  types  of  jet  planes. 
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Tabic  5.7.  Classification  Results:  A-10  Jet  No.  1  vs.  A-10  Jet  No.  2 


SNR 

Data 

Training 

Length 

A-10  No.  1 

No.  of  Correct 

Classifications 
(out  of  50) 

A-10  No.  2 

No.  of  Correct 

Classifications 
(out  of  50) 

Average  Time 
to  Classify 
A-10  No.  1 
(msec) 

Average  Time 
to  Classify 
A-10  No.  2 
(msec) 

40  dB 

6250 

49 

47 

10.01 

9.94 

2500 

11 

50 

9.68 

9.52 

20  dB 

6250 

43 

23 

9.68 

9.67 

2500 

19 

42 

9.65 

9.54 

10  dB 

6250 

43 

16 

9.54 

9.57 

2500 

15 

37 

9.52 

9.52 

This  is  obviously  a  more  difficult  problem  than  discriminating  between  a  jet  plane 
and  a  propellor  plane. 

When  comparing  the  different  jets,  the  results  were  more  balanced  than  those 
seen  comparing  the  jets  with  the  propellor  plane.  For  the  experiments  conducted 
comparing  Jet  No.  1  and  Jet  No.  2,  shown  in  Table  5.7,  percentages  of  correct  clas¬ 
sification  greater  than  eighty  percent  were  typical  for  about  half  of  the  scenarios. 
The  rest  of  the  results  ranged  from  0%  to  about  50%.  The  time  needed  for  the  al¬ 
gorithm  to  make  a  classification  decision  is  somewhat  higher  than  the  time  needed 
when  comparing  a  jet  plane  to  the  propellor  plane. 

Table  5.8  contains  the  results  for  the  experiments  conducted  comparing  the 
A-10  Jet  No.  1  with  Jet  No.  3.  From  this  table  one  can  see  that  the  algorithm  was 
able  to  successfully  classify  Jet  No.  1  correctly  for  the  environments  with  SNRs  of  20 
dB  and  10  dB,  while  Jet  No.  3  performed  poorly  under  these  conditions.  The  results 
were  opposite  for  the  case  of  an  SNR  of  40  dB.  Classifications  in  the  20  dB  and  10  dB 
cases  were  made  almost  immediately  after  reaching  the  minimum  allowable  length. 
For  the  40  dB  case  the  classification  time  was  significantly  longer. 

Results  for  comparing  the  A-10  Jet  No.  2  and  the  A-10  Jet  No.  3  are  contained 
in  Table  5.9.  From  this  table  one  can  see  that  nearly  perfect  results  for  classifying 
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Tabic  5.8.  Classification  Results:  A-10  Jet  No.  1  vs.  A-10  Jet  No.  3 


SNR 

Data 

Training 

Length 

A-10  No.  1 

No.  of  Correct 

Classifications 
(out  of  50) 

A-10  No.  3 

No.  of  Correct 

Classifications 
(out  of  50) 

Average  Time 
to  Classify 
A-10  No.  1 
(msec) 

Average  Time 
to  Classify 
A-10  No.  3 
(msec) 

40  dB 

6250 

11 

50 

10.56 

9.52 

2500 

5 

50 

9.53 

9.52 

20  dB 

6250 

50 

1 

9.52 

9.66 

2500 

50 

2 

9.52 

9.52 

10  dB 

6250 

50 

0 

9.52 

9.52 

2500 

50 

0 

9.52 

9.52 

Table  5.9. 


Classification  Results:  A-10  Jet  No.  2  vs. 


A-10  Jet  No.  3 


SNR 

Data 

Training 

Length 

A-10  No.  2 

No.  of  Correct 

Classifications 
(out  of  50) 

A-10  No.  3 

No.  of  Correct 

Classifications 
(out  of  50) 

Average  Time 
to  Classify 
A-10  No.  2 
(msec) 

Average  Time 
to  Classify 
A-10  No.  3 
(msec) 

40  dB 

6250 

40 

49 

12.48 

11.40 

2500 

50 

49 

13.39 

11.65 

20  dB 

6250 

49 

4 

9.53 

9.56 

2500 

50 

0 

9.52 

9.52 

10  dB 

6250 

50 

0 

9.52 

9.52 

2500 

50 

0 

9.52 

9.52 
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Jet  No.  2  were  obtained.  Classifications  occurred  almost  immediately  after  reaching 
full  order  for  the  20  dB  and  10  dB  cases.  While  more  time  was  needed  for  the  40 
dB  cases.  For  Jet  No.  3  nearly  perfect  classification  results  were  obtained  for  the  40 
dB  cases.  This  performance  fell  to  nearly  0%  for  the  20  dB  ad  10  dB  cases.  Similar 
to  the  experiments  with  the  propellor  plane  and  Jet  No.  1,  as  the  SNR  is  dropped 
significantly,  degradation  only  in  Jet  No.  3  tends  to  show  that  the  classifier  is  skewed 
towards  Jet  No.  2  questioning  whether  the  results  for  Jet  No.  2  at  SNRs  of  20  dB 
and  10  dB  are  reliable. 

In  another  set  of  experiments,  the  classifier  designed  using  the  propellor  and 
Jet  No.  1  was  tested  on  data  from  all  four  classes  (propellor,  Jet  No.  1,  Jet  No.  2  and 
Jet  No.  3).  When  tested  on  Jet  No.  2  and  Jet  No.  3,  the  classification  was  considered 
to  be  correct  if  the  classifier  chose  the  “Jet”  category,  and  incorrect  if  the  classifier 
classified  these  jet  planes  as  propellor  planes. 

Results  of  these  experiments  are  listed  in  Table  5.10.  From  these  results,  it 
can  be  seen  that  the  classifier  was  able  to  correctly  classify  the  A-10  Jets  No.  2  and 
No.  3  as  jets  perfectly  for  the  40  dB  scenarios.  When  the  SNR  was  reduced  to  20 
dB,  the  A-10  Jet  No.  3  was  classified  nearly  perfectly,  while  the  A-10  Jet  No.  2  had 
at  most  a  classification  rate  of  about  60%.  Under  the  10  dB  conditions,  the  classifier 
performed  poorly  when  trying  to  classify  the  A-10  Jets  No.  2  and  No.  3. 

2.  Analysis  of  Single-channel  vs.  Multirate  vs.  Multichannel 
Simulations 

Experiments  identical  to  those  conducted  for  the  multirate  data  were  con¬ 
ducted  for  the  single  channel  and  multichannel  scenarios.  In  most  situations,  it  was 
possible  to  adequately  classify  the  target  in  these  other  scenarios.  Therefore,  to  com¬ 
pare  the  effectiveness  of  the  methods,  the  average  number  of  data  samples  needed  to 
make  a  classification  were  also  compared.  Tables  5.11  through  5.13  show  results  for 
some  of  the  experiments  for  an  environment  with  an  SNR  of  40  dB.  Table  5.11  shows 
that  under  these  conditions,  the  classifier  was  able  to  correctly  classify  the  propellor 


110 


Table  5.10.  Classification  Results:  Class  0  =  Propellor  and  Class  1  =  A-10  Jet  No.  1 


SNR 

Data 

Training 

Length 

Propellor 

No.  of  Correct 

Classifications 
(out  of  50) 

Jet  No.  1 

No.  of  Correct 

Classifications 
(out  of  50) 

Jet  No.  2 

No.  of  Correct 

Classifications 
(out  of  50) 

Jet  No.  3 

No.  of  Correct 

Classifications 
(out  of  50) 

Average  Time 
to  Classify 
Propellor 
(msec) 

Average  Time 
to  Classify 

Jet  No.  1 
(msec) 

Average  Time 
to  Classify 

Jet  No.  2 
(msec) 

Average  Time 
to  Classify 

Jet  No.  3 
(msec) 

40  dB 

6250 

50 

50 

50 

50 

9.52 

9.52 

9.52 

9.52 

2500 

50 

50 

30 

50 

9.52 

9.52 

10.24 

9.52 

20  dB 

6250 

50 

50 

32 

47 

9.52 

9.52 

9.57 

9.53 

2500 

50 

39 

3 

41 

9.52 

9.75 

9.55 

9.54 

10  dB 

6250 

50 

21 

2 

9 

9.52 

9.52 

9.52 

9.52 

2500 

50 

10 

0 

10 

9.52 

9.52 

9.52 

9.53 

and  Jet  No.  2  in  all  these  scenarios.  In  addition,  it  can  be  seen  that  the  time  needed 
to  make  the  classification  is  roughly  the  same  for  the  single-channel,  multichannel 
and  multirate  cases. 

Table  5.11.  Classification  Results  for  Propcllor  and  Jet  No.  2 


Scenario 

Data 

Training 

Length 

Propellor 

No.  of  Correct 

Classifications 
(out  of  50) 

Jet  No.  2 

No.  of  Correct 

Classifications 
(out  of  50) 

Average  Time 
to  Classify 
Propellor 
(msec) 

Average  Time 
to  Classify 

Jet  No.  2 
(msec) 

Single-Channel 

6250 

50 

50 

10.12 

9.56 

2500 

50 

50 

9.56 

9.52 

Multichannel 

6250 

50 

50 

9.52 

9.52 

2500 

50 

50 

9.76 

10.28 

Multirate 

6250 

50 

50 

9.53 

9.52 

2500 

50 

50 

9.52 

9.53 

In  table  5.12,  it  can  be  seen  that  ability  to  distinguish  between  Jet  No.  1  and 
Jet  No.  2  becomes  more  difficult  than  between  the  propellor  plane  and  jets  for  all 
three  scenarios.  The  multirate  classifier,  however,  is  able  to  make  its  classifications 
in  a  shorter  amount  of  time. 

In  the  example  of  the  A-10  Jet  No.  2  versus  the  A-10  Jet  No.  3,  Table  5.13, 
the  classifiers  again  performed  fairly  well.  Although,  the  amount  of  time  necessary 
to  make  classification  decisions  for  all  three  cases  was  slightly  higher  than  for  some 
of  the  other  simulations,  the  multirate  classifier  on  average  needed  less  time  to  make 
its  classification. 

From  these  tables  it  can  be  seen  that  the  multirate  case  has  some  advantage 
over  the  single-channel  case.  In  low  noise  environments  where  the  single-channel, 
multirate  and  multichannel  scenarios  all  had  a  high  success  rate  of  proper  classifica¬ 
tion,  the  multirate  scenario  was  able  to  classify  the  proper  signal  in  a  shorter  amount 
of  time.  It  was  noticed  that  the  multichannel  channel  case  at  times  needed  more  time 
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Table  5.12.  Classification  Results  for  Jet  No.  1  and  Jet  No.  2 


Data 

Jet  No.  1 

Jet  No.  2 

Average  Time 

Average  Time 

Scenario 

Training 

No.  of  Correct 

No.  of  Correct 

to  Classify 

to  Classify 

Length 

Classifications 

Classifications 

Jet  No.  1 

Jet  No.  2 

(out  of  50) 

(out  of  50) 

(msec) 

(msec) 

Single-Channel 

6250 

50 

36 

9.90 

35.78 

2500 

45 

49 

15.96 

22.78 

Multichannel 

6250 

50 

35 

9.89 

32.42 

2500 

46 

50 

15.37 

22.51 

Multirate 

6250 

49 

47 

10.02 

9.94 

2500 

11 

50 

9.68 

9.52 

Table  5.13.  Classification  Results  for  Jet  No.  2  and  Jet  No.  3 


Data 

Jet  No.  2 

Jet  No.  3 

Average  Time 

Average  Time 

Scenario 

Training 

No.  of  Correct 

No.  of  Correct 

to  Classify 

to  Classify 

Length 

Classifications 

Classifications 

Jet  No.  2 

Jet  No.  3 

(out  of  50) 

(out  of  50) 

(msec) 

(msec) 

Single-Channel 

6250 

35 

50 

20.10 

10.74 

2500 

50 

49 

15.13 

13.08 

Multichannel 

6250 

37 

50 

17.29 

11.31 

2500 

50 

48 

12.05 

15.18 

Multirate 

6250 

40 

49 

12.48 

11.40 

2500 

50 

49 

13.48 

11.65 
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to  make  a  classification  decision  than  the  multirate  case.  Although  there  is  more 
data  available  in  the  multichannel  scenario,  the  offset  observations  of  the  multirate 
scenario  adds  more  statistical  information  to  the  solution. 

3.  Key  Parameters  for  the  Multirate  Case 

There  were  three  main  characteristics  of  the  data  that  were  analyzed  for  their 
effect  on  the  classification  performance.  They  were  the  dominant  frequency  of  the 
target  data,  the  signal-to-noise  ratio  and  the  training  data  length.  In  addition,  the 
performance  between  single-channel,  multirate  and  multichannel  observations  were 
compared.  Each  is  described  individually  in  the  following  subsubsections,  however, 
their  impact  is  not  independent  of  each  other.  So,  if  noticeable,  their  relative  impact 
on  each  other  is  discussed. 

a.  Dominant  Frequency 

While  analyzing  the  data,  it  became  apparent  that  the  classifier  per¬ 
formed  much  better  when  one  target  data  class  had  a  dominant  frequency  higher  than 
that  of  the  other  class  to  which  it  was  being  compared.  This  was  especially  noticeable 
in  the  high  noise  environment  scenarios.  Under  the  high  noise  conditions,  the  jets 
were  not  classified  correctly  with  any  degree  of  success  against  the  propcllor  plane. 
Only  under  the  low  noise  condition  ( SNR  =  10  dB)  did  the  jets  display  moderate 
success  with  classification  against  the  propellor  plane.  Likewise,  Jet  No.  1  had  a 
higher  correct  classification  rate  when  attempting  to  classify  against  Jets  No.  2  and 
No.  3.  And,  Jet  No.  2  performed  better  against  Jet  No.  3. 

b.  Signal-to- Noise  Ratio 

By  analyzing  the  simulation  data,  it  became  noticeable  that  under  low 
noise  conditions  ( SNR  =  40  dB),  the  three  A-10  jets  were  classified  properly  against 
the  propcllor  plane.  As  the  training  sample  length  was  reduced,  the  performance 
degraded  dramatically.  In  addition,  the  higher  noise  conditions  reduced  the 
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performance  of  the  jets  significantly  when  compared  to  the  propcllor  plane,  ft  appears 
that  the  dominant  frequency  has  a  large  effect  on  the  outcome  of  the  classifier,  in  any 
condition  other  than  ideal. 

c.  Training  Data  Length 

As  might  be  anticipated,  the  longer  the  training  data  length,  the  better 
the  classifier  typically  performed.  Not  only  did  the  classification  rate  generally  in¬ 
crease,  but  the  number  of  samples  needed  to  make  a  proper  classification  was  reduced. 
For  the  audio  data  used,  however,  when  the  size  of  the  training  data  was  increased  to 
12500  samples  the  performance  of  the  classifier  actually  degraded.  It  was  determined 
that  at  this  length  the  data  available  for  testing  became  statistically  different  enough 
from  the  training  data  and  caused  a  degradation  in  the  performance  of  the  classifiers. 

G.  CONCLUSION 

In  this  chapter,  it  was  shown  that  with  some  modifications  to  the  single- 
channel  classifier,  one  can  implement  a  sequential  classifier  that  accepts  multirate 
data.  In  addition,  using  the  multirate  Levinson  recursion  developed  in  Chapter  IV; 
one  can  derive  the  needed  class  parameters  for  the  multirate  sequential  classifier.  The 
multirate  sequential  classifier  can  handle  multiple  channels  of  observations  occurring 
at  different  sampling  intervals. 

Simulations  were  conducted  to  analyze  the  effectiveness  of  the  multirate  se¬ 
quential  classifier  compared  with  a  single-channel  classifier  and  a  multichannel  clas¬ 
sifier  where  all  channels  were  sampled  at  the  same  rate.  It  was  seen  that  under  low 
noise  environments  all  three  classifiers  performed  well,  but  the  multirate  classifier  was 
able  to  make  classification  decisions  in  a  shorter  amount  of  time  due  to  the  added 
information  from  the  offset  observations.  The  amount  of  addiitonal  information  con¬ 
tained  in  the  multirate  classifier  is  dependent  upon  the  disparity  between  channel 
sampling  rates. 
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VI.  CONCLUSIONS  AND  FUTURE  WORK 


The  goal  of  this  work  was  to  further  develop  the  theory  and  applications  of 
statistical  multirate  signal  processing.  A  summary  of  the  work  and  contributions  of 
this  thesis  is  provided  here.  This  is  followed  by  a  section  of  suggestions  for  further 
research. 

A.  CONCLUSIONS 

In  Chapter  II,  the  foundation  necessary  to  describe  multirate  systems  is  es¬ 
tablished.  Part  of  this  foundation  requires  defining  some  key  terms  applicable  to 
multirate  systems,  such  as  fundamental  rate,  system  period  and  system  phase.  These 
terms  are  used  to  explicitly  describe  the  periodic  nature  of  the  multirate  system 
and  allow  for  the  appropriate  selection  of  values  for  the  periodic  components  in  the 
multirate  system.  The  definitions  for  describing  the  statistical  characterizations  of 
multirate  signals  and  building  blocks  (e.g.,  decimator  and  expander)  are  provided. 
Compact  matrix  forms  of  the  decimator,  expander  and  filters,  which  take  advantage 
of  Kronecker  product  relations,  are  presented.  Finally,  the  LPTV  filter  bank  of  [Ref. 
17]  is  generalized  to  apply  to  systems  with  different  input  and  output  rates. 

In  Chapter  III,  the  explicit  direct  form  of  the  multirate  optimal  estimator  is 
developed.  This  optimal  estimator  uses  multiple-input  signals  observed  at  different 
sampling  rates  to  estimate  a  desired  signal.  In  addition,  a  recursive  innovations 
form  of  this  multirate  optimal  estimator  is  derived.  This  innovations  form  of  the 
optimal  estimator  separates  the  direct  form  optimal  filters  into  modified  optimal 
filters  and  cross  filters.  These  cross  filters  remove  any  information  from  one  signal 
that  is  contained  in  the  other  signals,  in  an  ordered  fashion.  In  essence,  a  new  set 
of  input  signals  that  are  mutually  orthogonal  are  derived  and  used  as  the  inputs  to 
the  modified  optimal  filters.  Using  the  recursive  form  of  the  optimal  filter  allows 
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one  to  calculate  the  relative  change  in  performance  in  the  optimal  estimator  when 
input  signals  are  added  or  removed  from  the  system.  Through  simulations,  it  is 
demonstrated  that  optimal  filtering  of  multiple  multirate  channels  for  an  AR  process 
and  a  sinusoidal  process  can  provide  improved  performance  over  optimal  filtering 
using  a  single  channel,  even  if  the  secondary  channel  when  used  by  itself  has  high 
error  variances. 

In  Chapter  IV,  the  optimal  filtering  problem  is  specialized  for  the  case  of 
optimal  multirate  linear  prediction.  The  multirate  Normal  equations  are  derived 
for  a  multirate  system  with  multiple  input  and  output  signals,  which  are  observed  at 
different  sampling  rates.  For  a  multirate  system  with  a  system  period  K,  there  are  up 
to  K  distinct  sets  of  prediction  coefficients  and  error  covariance  matrices  that  apply  in 
a  periodic  fashion.  An  efficient  method  for  calculating  the  multirate  linear  prediction 
coefficients  and  error  variances  is  developed  through  the  use  of  the  multichannel 
Levinson  recursion  and  generalized  triangular  UL  factorization.  A  specific  algorithm 
for  the  generalized  triangular  factorization  is  included  in  Appendix  E. 

In  Chapter  V,  a  multirate  sequential  classifier  is  derived  starting  from  the 
basic  theory  of  sequential  hypothesis  testing.  It  is  shown  that  classifier  parameters 
needed  for  implementing  the  multirate  sequential  classifier  are  the  same  as  those  for 
multirate  linear  prediction.  A  multirate  sequential  classifier  is  then  implemented  and 
tested  using  audio  hies  of  a  propcllor  plane  and  three  A-10  jet  aircraft.  The  experi¬ 
ments  tested  the  classifier  performance  in  selecting  between  the  propellor  plane  and 
jet  aircraft  as  certain  system  parameters  are  changed.  These  parameters  consisted  of 
the  signal-to-noise  ratios  of  the  observed  signal  and  the  length  of  the  training  data. 
In  addition,  the  performance  of  the  multirate  classifier  is  compared  to  that  of  the 
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single-channel  and  multichannel  classifiers  using  similar  data.  These  experiments 
showed  that  improved  performance  can  be  obtained  by  using  two  channels  of  obser¬ 
vation  data  even  if  one  channel  is  sampled  at  a  lower  rate. 

B.  TOPICS  FOR  FURTHER  INVESTIGATION 

Multirate  signal  processing  continues  to  be  an  active  area  of  research  and  many 
opportunities  exist  for  further  research,  particularly  in  the  area  of  statistical  signal 
processing.  The  optimal  multirate  filter  developed  in  this  dissertation  is  an  FIR 
filter.  The  infinite  impulse  response  (HR)  filter  was  not  considered  here;  however, 
multirate  theory  applies  equally  well  to  the  HR  case  and  also  would  provide  significant 
contributions  to  the  theory  of  multirate  signal  processing. 

Another  opportunity  for  further  work  on  multirate  systems  is  in  the  area  of 
lattice  structures.  Some  work  on  lattice  structures  has  already  been  researched,  but 
extensive  development  of  lattice  theory  as  it  applies  to  multirate  linear  prediction  can 
lead  to  more  convenient  representation  and  implementation.  In  addition,  this  disser¬ 
tation  focused  on  optimal  solutions  using  time-domain  characteristics.  Solutions  to 
optimal  filtering  of  multirate  systems  using  frequency-domain  techniques  would  fur¬ 
ther  expand  the  theory  and  utility  of  multirate  statistical  signal  processing.  Research 
into  filter  design  of  periodically  correlated  scalar  time  series  using  a  spectral  factor¬ 
ization  algorithm  for  wide-sense  stationary  vector  processes  has  already  been  started 
by  Spurbeck  and  Scharf  [Ref.  29] . 

Research  from  a  companion  thesis  [Ref.  50]  has  extended  the  concepts  of  opti¬ 
mal  multirate  filtering  to  two-dimensional  (2-D)  signal  processing  and  high-resolution 
image  reconstruction.  The  reader  of  this  thesis  may  want  to  refer  to  [Ref.  50]  for 
some  additional  topics  not  treated  here. 
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APPENDIX  A.  POWER  SPECTRAL  DENSITY 
FOR  MULTIRATE  PROCESSES 


This  appendix  contains  the  derivations  of  the  power  spectral  density  functions 
for  the  the  decimator,  expander  and  filters  presented  in  Chapter  II.  The  standard 
equation  for  the  power  spectral  density  function  of  a  stationary  process 

OO 

SM  =  X  (A.i) 

/=— OO 

can  be  applied  to  the  input  and  output  correlations  independently.  However,  any 
difference  between  sampling  rates  of  the  input  and  output  signals  makes  it  impossible 
to  develop  a  power  spectral  density  function  for  the  cross-correlation  dependent  on 
a  lag  only.  Therefore  it  is  necessary  to  use  a  two-dimensional  power  spectral  density 
function.  This  function  is  defined  as 

OO  OO 

S'2D(A1,  A0)  =  J5  51  R[n\,  n0]e-jAinie-jA°n°.  (A.2) 

m=— oo  riQ= — oo 

For  the  input  to  a  decimator,  expander  or  filter,  the  power  spectral  density  is 
defined  as 

OO  OO 

S2XD{ A!,Ao)=  55  55  Rx[ni,n0]e-^e-^n°.  (A.3) 

Til  — — oo  TlQ  = — oo 

Using  the  time-lag  form  (when  appropriate),  the  power  spectral  density  function 
becomes 

OO  OO 

SlD(\,u)=  Y,  Y  (A. 4) 

n=— oo  l=— oo 

which  reduces  to  (A.l)  when  the  signal  is  wide  sense  stationary. 
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A.  DECIMATION 


The  operation  of  decimation  is  depicted  in  Fig.  A.l  and  defined  as 


y[my]  =  x[Lmy\  for  m  —  . . . ,  —  1,  0, 1, . 


'K] 


II 


’[my] 


Figure  A.l.  Decimator 

The  correlation  functions  for  the  decimator  are  as  shown  in  Table  2.2  of  Chap¬ 
ter  II.  The  output  power  spectral  density  function  for  decimation  is 


OO  OO 


5f(A,,A„)=  V  V  Ry{m9,rn'y]e-»'m>e-»°m'- 


my=—oc  m  y  = — oo 
OO  oo 


E  E  Rx[Lmy,  Lm'y 


i'y^- 


OO  OO 


^2  Rx[LmyiLm'y\e~jy/'ll"^'"'ye~J^l±J)1J"Ly 


e~j\im'y  e-j\0my 


,-j(*l/L)Lm'  -j(\0/L)Lmy 


mv=— oo  m'=— oo 


n2Df\  \  \  _  q2D  I  ^1  A0 

by  fAi,  A0)  -  ax 


In  the  time-lag  representation,  the  output  power  spectral  density  function  is  defined 
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as 


oo  oo 

Sf( A,«)=  X]  E 

my— — oo  ly=— oo 
oo  oo 

=  ]T  Rx[L'my]  Lly\e~^Xrny e~^lv 

TTly —  OO  l  y  = — OO 

OO  oo 

my= — OO  ly=— oo 


5f(A,a,)=5r(p^ 


or  in  lag  representation  only 


OO 

S„M  =  E  ^[ye"*"' 

ly  —  OO 

oo 

=  ]T 

ly -  OO 

OO 

=  ^  iqLye-j(aj/L)i^ 

ly  —  OO 

S„M  =  5.  Q 


The  cross-power  spectral  density  function  is  given  by 

OO  OO 

S™(  \x,\y)=  J2  Y  Rxy[mx,my\e-ix°m°e-iXymy 

mx=— oo  my=— oo 
oo  oo 

=  Y  Y  Rx[mx,Imiv]e-j^e-j^ 

mx= — oo  my=— oo 
oo  oo 

=  ^  -Rr  [mx,  Lmy\e^Xxrnxe~^Xy^L^Lrny 

mx=— oo  vriy— — oo 
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In  the  time-lag  representation,  the  cross-power  spectral  density  function  is  defined  as 

OO  OO 

S™(\x:u,)=  Y.  'll  RXy{mx;lv}e-ix-m-e-^‘- 

mx=—oo  ly=— oo 
OO  oo 

=  J2  Rx[mx;Lly-{L-l)mx]e-jX^e-j^ 

mx —  oo  l  y  —  oo 


oo  oo 

^  ^  ^  ^  Rx\j^Xl  Lly 

TYly —  OO  l  y  —  OO 

OO  oo 

^  ^  ^  ^  ^  r  J^  X  ■  Lly 

my=— oo  ly=—oo 


S™(\X,Uy)  =  S, 


xy 


2D 

x 


L  1, ,  LJy\ 

—Uv'T) 


_ \)mx]e~jXxmx e  ^(L£la;!/)ma:£— i(  LL1^y)mx 


B.  EXPANSION 


The  operation  of  expansion  is  depicted  in  Fig.  A. 2  and  defined  as 


y[m\  = 


x[m/I ]  m  div  I 


otherwise 


AmA 


t  / 


Figure  A. 2.  Expander 


The  correlation  functions  for  the  expander  are  as  shown  in  Table  2.3  in  Chapter 
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II.  The  output  power  spectral  density  function  for  expansion  is 


OO  OO 


SyD(Xi,  A0)  =  ^  ^  Ry[my,  m'y\e-j^mve-jXomy 


mv=— oo  m'=— oo 


E  E  Rx[my/I,m'y/I}e-^mye-^my 


mv=— oo  m'=— oo 


oo  oo 


mv=—oo  m'=— oo 


^2D(Ai,A0)  =  S'2d  (/Ai,  JA0) 


In  the  time-lag  representation,  the  output  power  spectral  density  function  is  defined 


SyD( \,u)=  Y  Ry[myilv 


e  “e 


TTl/ij —  OO  ly —  OO 


E  E 


TTXij —  OO  l  y —  OO 


E  E  Rx[my/I'i  ly/I]1 


-j(I\)my/Ie-j(Iu)ly/I 


SlD{\u)=SlD  (/A,  Iu) 

A  lag  representation  for  the  expander  does  not  exist,  since  the  output  is  not  wide 
sense  stationary.  The  cross-power  spectral  density  function  is  given  by 


OO  OO 


S%(\x,  Xy)=  Y  E  RxV[mx,mv]e-j^e-j^ 


mx — — oo  mv=— oo 


Y  Y  Rx[mx,mv/I]e-jX*m*e-jX*m* 


mx — — oo  mv=— oo 


oo  oo 


y  E  jAxmie 


-j\xmx—j(I\y)m.y/I 


mx=— oo  mv=—oo 


S2°(\x,\y)  =  S2xD[\xY 
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In  the  time-lag  representation,  the  cross-power  spectral  density  function  is  defined  as 

OO  OO 

S™(\x:u,)=  Y.  'll  Rx,{mx;lv}e-ix-m-e-^‘- 

mx=—oo  ly=— oo 
OO  oo 

=  J2  E  Rx[mx-(ly  +  (I-l)mx)/I]e-jX^e-j^ 

mx —  oo  l  y  —  oo 
oo  oo 

Rx[mx ;  (ly  +  (/  -  l)mx)/I]e-R^e-K-{I-1)uJymx)e-KIuJy)ly/Ie-KI-1)uJymx 

ITT. *y —  OO  l  y  —  OO 

OO  oo 

Rx[mx-  (ly  +  (/  - 

my=— OO  ly=— oo 

=  SX2D  (A.  -  (/  -  lK,/^) 


C.  FILTERS 

Filtering  is  depicted  in  Fig.  A. 3  and  defined  as 

OO 

y[m]  —  h[r\x[m  —  r]. 


x\m\ 


y[m\ 


Figure  A. 3.  Filter 


The  correlation  functions  for  the  filter  are  as  shown  in  Table  2.4  in  Chapter 
II.  The  output  power  spectral  density  function  for  a  filter  is 

OO  OO 

S2yD( Al5Ao)=  ^  ^  Ry[m,rn'}e-jX'me-jXom' 

m=— oo  m'=—oo 
oo  oo 

=  m']  *  /i[m]  *  h*[m'])  e~^Xime^Xorn' 

m=— oo  m'=— oo 

Af  (A,,A0)  =  5^(Ai,Ao)^(A1)^*(A0) 
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In  the  time-lag  representation,  the  output  power  spectral  density  function  is  defined 


as 

OO  OO 

sf( A,w)  =  X 

m=—oo  l=— oo 
OO  oo 

=  ^  ^  /]  *  h[m]  *  /?*[— /])  e-^Xme-^1 

m= — oo  l=—oc 

S2D(\,u;)  =  S2xD(\,u)H(\)H*(u) 

or  in  lag  representation  only 

OO 

s»  =  ]r 

l=— OO 
oo 

=  ]T  (i?x[Z]*/l[Z]*/l*[Z])e-^ 

/=— OO 

=  Sx{lu)H(u;)H*(u;) 

The  cross-power  spectral  density  function  is  given  by 

OO  OO 

O  V,Ao)=  V  v 

m=— oo  m'=— oo 
oo  oo 

m=— oo  m'=— oo 

S^(Ai,Ao)  =  S“(Ai,\,)/f*(Ao) 

In  the  time-lag  representation,  the  cross-power  spectral  density  function  is  defined  as 

OO  OO 

S“(A.V)=  V  £ 

m— — oo  /=— oo 
oo  oo 

=  X]  i2x[m;Z]  * /i*[-Z]e-J'Ame“J'wZ 

m=— oo  Z=— oo 

S^(A,u;)  =  ^(A,u;)fr(u;) 
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D.  SUMMARY  OF  CORRELATION  AND  POWER  SPEC¬ 
TRAL  DENSITY 

For  ease  of  reference,  a  summary  of  the  correlation  functions  and  power  spec¬ 
tral  densities  for  decimation,  expansion  and  filtering  is  provided  in  Tables  A.l  through 
A. 3. 


Table  A.l.  Summary  of  Decimation 


:  [  mx  ]  - 


iL 


*y[my] 


Ry[my,m'y]  =  Rx[Lmy,  Lm'y\ 

Sr(A1,A0)  =  SjD(^,^) 

Ry\lTly^  /y]  RX  \LlTly  ,  LI  y] 

Sf(  \,u)  =  S™  (A,  f) 

Ry[ly]  =  Rx[Lly\ 

SM  =  S,  (|) 

RXy[mx,  my\  =  Rx[mx,  Lmy] 

=  (A 

Rxy\jYlx’)  ly\  Rx\jYlxi  Lly  i^L  1^777/^] 

Sxy  (Ax,  ujy)  —  Sx  ^\x  +  ^  UJyi  L  ) 
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Table  A. 2.  Summary  of  Expansion 


*K] 


t  / 


S3 

C£0 

II 

Rx[my/I ,  m'y/I\  my,  m'y  div  I 

0  otherwise 

^2D(A1,Ao)  =  S’^(/A1,/A0) 

S3 

II 

Rx[mv/I;ly/I \  my,  ly  div  I 
]  otherwise 

SlD(\u)  =  SlD(I\Iu>) 

Rxy  •>  ^'y] 

j  Rx[mx,  my/I\  my  div  I 

1  0  otherwise 

S™(\x,\y)  =  SlD(\xJ\y) 

S3 

« 

W 

II 

Rx[mx]  (ly  +  (I  -  1  )mx)/I]  mx  -  ly  div  / 

0  otherwise 

S™(K,coy)  =  S2XD  (\x  -(I-  IK,  IiUy) 

Table  A. 3.  Summary  of  Filtering 
x[m\ - »  h[m]  - *  y[m] 


Ry[m,  m']  =  Rx[m,  m']  *  h[m]  *  h*[m'] 

S2yD( Ar,A0)  =  SxD(\i,  \0)H(\i)H*(\0) 

Ry[l]  =  Ry[l]*h[l]*h*[-l] 

Sy(u)  =  Sx(lu)H(u;)H*(u;) 

Rxy[m ,  m']  =  Rx[m,  vnl\  *  h*[m'] 

S™(\1,\o)  =  S2xd(\1,\0)H*(\0) 

Rxy[m,  1}  =  Rx[m ;  l]  *  h*[-l\ 

S™(\u\„)  =  S2n(\,v)H‘(u,) 
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APPENDIX  B.  KRONECKER  PRODUCT  AND 
REVERSAL  NOTATION 


A.  KRONECKER  PRODUCT 

This  appendix  provides  the  definition  of  a  Kronecker  Product  and  lists  some 
useful  relations.  A  thorough  discussion  of  Kronecker  products  and  matrix  calculus  is 
contained  in  [Ref.  62]  and  a  shorter  concise  synopsis  is  contained  in  [Ref.  63]. 

Given  two  matrices  A  and  B.  the  Kronecker  product  A  ®  B  is  defined  as 

anB  012B  . . .  aijB 

a2iB  CI22B  . . .  a2,B 

A  = 

aaB  a,;2B  . . .  a^B 

Some  basic  properties  and  relations  for  Kronecker  products  are  provided  in  Table  B.l. 

Table  B.l.  Some  Properties  of  Kronecker  Products 

(1)  (A  ®  B)  (D  <g)  G)  =  AD  <g)  BG 

(2)  (A  +  H}<g)(B  +  R)  =  A<g)B  +  A<g)R  +  H<g>B  +  H<g)R 

(3)  (A®B)t  =  At®Bt 

(4)  (A  <g>  B)_1  =  A-1  <g>  B”1 

(5)  tr  (A  <g)  B)  =  tr  (A)  •  tr  (B) 

(6)  (A  <g>  B)  «g)  C  =  A  ®  (B  <g>  C) 
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B.  REVERSAL  NOTATION 


Given  a  vector  v  of  length  P  whose  elements  are 

r  i  T 


a  i  CL'2 


Up 


the  reversal  of  the  vector,  denoted  by  v,  is  defined  by 

r  iT 


v  = 


U  p 


Up- 1 


CL\ 


Given  a  P  x  Q  matrix  A  with  elements 


flip 

Ol,2  ‘  ‘  ' 

■  Oi,Q 

A 

= 

02,1 

02,2  '  '  ' 

■  o2,q 

Op,  1 

Op, 2  '  '  ' 

'  Op,Q 

the  reversal  of  A  is  given  by 

[ 

aP,Q 

OPQ-1 

r  aP-l,Q  1,Q  — 1  •  •  •  dp- 1,1 

— 

al,Q  al,Q-l  '  '  '  °1,1 

A  few  important  properties  of  vector  and  matrix  reversal  are  listed  in  Table  B.2. 
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Table  B.2.  Properties  of  Reversal  (from  [Ref.  18]) 


Quantity 

Reversal 

Matrix  product 

AB 

AB 

Matrix  inverse 

A"1 

(A)-1 

Matrix  conjugate 

A* 

(A)* 

Matrix  transpose 

AT 

(A)T 
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APPENDIX  C.  DERIVATION  OF 
PARAMETERS  FOR  INNOVATIONS  FORM  OF 
THE  OPTIMAL  FILTER 


In  Chapter  III  recursive  optimal  filters  were  used  to  present  the  optimal  filter 
in  its  innovations  representation.  This  appendix  contains  the  derivations  related  to 
these  optimal  filters. 


A.  REPRESENTING  THE  OPTIMAL  FILTER  IN  A  RE¬ 
CURSIVE  FORM 


The  Wiener-Hopf  equations  pertaining  to  an  optimal  multirate  filter  with  M 
input  channels  are  given  by  (3.7).  The  solution  for  the  filter  coefficients  can  be  written 


as 


where 


and 


bM 


lb(fe) 

cZ(M+l)  ’ 


p  (fc) 

'*11 

p  C') 
^lM 

p  (k) 

b(fc) 

a(M+l) 


In  general,  for  the  ith  signal  (C.2)  and  (C.3)  can  be  defined  as 


(C.l) 


(C.2) 


(C.3) 


(C.4) 


R 


(fc)  *T 
'il 


R 


(fc) 
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and 


b(fc)  - 

udi  ~ 


s(fc) 

■dl 


;(fc) 

d(2  — 1) 


In  order  to  derive  a  recursive  form  of  the  optimal  filter,  the  inverse 
is  hrst  written  in  a  recursive  form.  By  partitioning  into 


c[k)  = 


/■l(k) 

'-'i-l 

■p>  (fc)  *T 

13  a 

R.W 

where 


B(fc)  = 

o 


R 


(fc) 

ij 


"R 


the  recursive  form  of  C (see  [Ref.  69])  can  be  expressed  as 


c(fc)_1  = 


1 

'~>i— 1 

1 

o 

+ 

1 - 

o 

0 

1 

^i- 1 

i 

o 

+ 

1 - 

o 

0 

-G 

I 


(fc) 


E 


(fc)- 


_ j 


(fe)lj,(fc)  1  *T  -G^E^  1 


G-  E-  G 

_ gW  1Q(k)*T 


i  i 


E 


(fc)- 


where 


Q(fc)  =  cw 


and 


Ef>  =  R|‘>  -  Bj*,rCg"6f 

_  dW  _  r(fc)*TD(fc) 


(C.5) 

matrix  of 

(C.6) 

(C.7) 

(C.8) 

(C.9) 
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Applying  the  result  of  (C.7)  specifcally  to  the  Mth  signal  by  substituting  (C.7) 
into  (C.l)  leads  to 


K 


,(*) 


UM- 1  U 

0  0 


1  r-\(k)*T 


p(k)pi(k)  1 


■pi  (fc)  1  p 

■*-JM  'JM 


Separating  the  two  components  of  C®  then  produces 

hSfcr 


(fc)  *T 


E 


(fc)- 

M 


[b(fc)l 

udM 

f(fc) 

1  dM 

,(*) 

*M 


p(fc)  1  n 

W/- 1  U 

0  0 


b(fc) 

udM 

f(fc) 

+ 

1dM_ 

(fc)-rp(fc)  1  n  (k)  *T  pWi 7i (fc)  1 


pWpiW  p 
^M  ^M  ^M 

_ pi(fc)  1f~\(k)*T 

rjM  ^  t 


GW  ITU 
M 


M 


E 


(fcy 

M 


[b(fc)l 

udM 

f(fc) 

1  dM_ 

Extracting  the  terms  associated  with  h;^ ,  the  filter  is  found  to  be  equiv¬ 


alent  to 


uW  _  _ pi(7')  1p(^)*^'u(^)  _l_  ■p'W  1p(^) 

11 M  —  uclM  '  ^  M  1  dM 


_  piW  1  ( _  pt(fc)*T,  (fc) 

—  “M  1  1  dM  udM 


(C.10) 


and  the  rest  of  the  matrices  reduce  to 


dfe) 

*M-1 


_ p(fc)  1p,(^)  I 

—  '~'M-1  udM  “T 


p(k)-pi(7’)  1p(k)*7’ 

^M  'JM 


p  Wpi(fe)  1 

"'JM  “M 


;(*) 

1  dM 

s(fc) 

dM 


or 


h 


(fc) 

M—l 


pW  KW  i  p 

W/-1  DdM  + 


pi  w  1  p 

M  ^M  ^ 


■<(&)  'EW  _i_  pWpW  1  p(^)  *^7|(k)  _ pwp 

udM 


(fe)TTl(fe)  1~(fc) 


M 


M  M 


LdM- 


Further  manipulation  leads  to 


h 


(fc) 

M—l 


p(fc)  1r(fc) 
^ M—l  udM 


p  (fe)-p(fe) 

W/  ^M 


-G 


(fc)*rr(fc)  \ 
M  udM  J  i 


(C.ll) 
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and  substituting  (C.8)  and  (C.10)  into  (C.ll)  produces 


(fe) 


hi 


Ak) 
1M- 1 


_ p(fe)  1C(k)  p-t— 1  to (fe)  T,(fc) 

—  ^M-l  udM  • 


(C.12) 


Equation  (C.12)  can  also  be  expressed  as 

hifc) 


h(fc) 

“m-i 


_  p(fe)  1 

~  1 


b(fc) 


p-i 

Vfcf-1 


rd(iW-l) 

and  substituting  (C.7)  into  (C.13)  results  in 
h<fc> 


■g(fc) 

■o(fc) 


h(fc) 


(C.13) 


h(fc) 

11 M- 1 


p(fc)  1  n 
bM- 2  U 

0  0 


+ 


p  (fe) 


rj"  r-l  ^M-l 


'M 

■p(fe)_1p 


_p(fe)  rp(fc)  1 


(fe)  *T 
M— 1 


-l 


E^') 


’b(fc) 

cZ(M— 1) 

frf(M-l) 

p(fc)  1  Q 

UM- 2  U 

0  0 


+ 


(~i(k)  tiW 


(~t(k)  pi(fe)  1 


Tp(fe)  ‘pw« 

1  ^M-l 


(fe)  *T 


Extracting  the  terms  associated  with  h^_, ,  leads  to 

■CjM-lVjrM-l  ud(M- 1) 


E 


-i 


(fe) 

M—l 


B 


(fe) 


,(fc) 


P>(fe)  uW 


h(fc)  - 

uAf— 1  — 


p(fe)  ^(fe) 

1 1  d(M—l) 


Tfik)-1  r(k)*T-n(k) 
d(M- 


1  (fe) 


pC')'  p 

^M-l  1 V 


(fe)  u(fc) 


This  equation  can  be  expressed  equivalently  as 

u(fe)  _  E(fc)_1  7 p 
nM- 1  —  ^M-l  lr 


;(fc) 

d(M—l) 


p(fc)*Tr(fc) 
^M-l  ud(M-l) 


_  TP  (fe)  1  /  TO  l 


.(fe) 


P  (fe)  (fe) 

^M-l  d(M-1)M  I  “M 


(fe) 


(C.14) 
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The  terms  for  through  h®_2 


are  then  written  as 


h 


(fc) 

M—2 


c 


(fc)_1b 

M—2  u 


(fc) 

d(M-l) 


+ 


/~1  W  pO) 
^M-l^M-l  ''JM— 1 


/-<(&)  ph^')  1 

^M-l^M-l 


:(fc) 

3d(M-l) 

xW 

a(M— 1) 


_ /~<W  1‘p>(k)  n 

^M-2  J->(M-1)M11M 


(fc) 


/■iW  xfi(^)  -p(fc)  1 

^M-l^M-l  ^M-l  ^M-l^M-l 


B 


(fc) 


,(fc) 


(M-1)M“M 

■£>(fc)  u(fc) 
A*'(M-1)M11M 


This  equation  can  be  rearranged,  producing 


h!fc) 

h(fc) 

nM-2 


c 


(*)_1b 

M—2 


(fc) 

d(M-l) 


Tp(^)  1 
^M-l-^M-l 


(^(M-l) 


-G 


(fc)*rr(fc) 

M— 1  ud(M-l) 


+  G 


(*) 
M— 1 


.W 

M- 


R 


(fc) 

(M-l)M 


W/-1  °(M-1)M 


-c 


(fc) 


-1 


M—2 


D(M-1)MU 


(fc) 

M  • 


(C.15) 


Now  substituting  (C.8)  and  (C.14)  into  (C.15)  yields 


h) 


(fc) 


h(k) 

M—2 


—  ^ M—2  ud(M- 1) 


‘dW  u(^) 

^M-2  ■d(M-1)(M-1)11M-1 


/~<(k)  uW 

'-'M-2  J->(M-1)M11M  • 


(C.16) 


By  examining  the  pattern  of  how  (C.l)  is  separated  into  (C.10)  and  (C.12)  and  how 
(C.12)  is  separated  into  (C.14)  and  (C.16),  the  following  recursions  can  be  observed 


E 


(fc)- 


E 


(fc)- 

M 


(f(fc) 

y-di 

-  G!‘>*Tb«) 

M 

E 

1 

^ _ " 

i 

a' 

i=»+i 

PdM 

n(k)*TC(k) 

udM 

j(fc)*Tg(fc)\  j^(fc) 


(1  <  i  <  M  -  1) 


(C.17) 

(C.18) 
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where  the  matrices  and  E|A)  are  given  by 


G(,k)  =  { 


and 


■p  (fc) 

T>  (fe) 

■*^12 

"R 

^10-1) 

p(fc)*T 

rv12 

5s 

^ '  CNI 

■o(fc) 

B(fc) 

ll 

^2(1-1)  '  ' 

1 

l 

gl 

Ejfc; 1  = 

|  R-i? 

i  = 

1 

-l 


i  =  1 


1  <i<  M 


R-?fe)  -  G-fc)  *tB^  1  <i<M 


The  hlter  in  (C.17)  and  (C.18)  can  be  further  simplified  by  defining 


(fc)  _  -p(fc)-1  G(fc)*Tg(fc) 


H£'  =  E, 
h' (fc)  =  E 


a(fc) 


V 


G 


(k)*TC(k) 
i  udi 


(C.19) 

(C.20) 


and  substituting  these  terms  into  (C.17)  and  (C.18) 

M 

hf}  =  h{{k)  -  Y  (1  <  i  <  M  -  1)  (C.21) 


j=i+l 

h  S’  =  h'W  (C.22) 

The  hlter  coefficients  of  (C.21)  for  an  M-signal  system  can  be  expressed  in 
terms  of  the  cross  terms  and  orthogonal  filters  only,  which  is  useful  for 
development  of  the  recursive  error  variance.  Beginning  with  the  general  equation  for 
the  filters  for  signals  1  to  M  —  1,  (C.21),  but  using  the  index  i\  instead  of  j  for  the 
summation 

M 

h^  =  ht{k)+  Y  (-Hg^hf  (1  <  i  <  M  —  1)  (C.23) 

h=i+ 1 


the  hlter  for  signal  i  can  be  expanded  by  substituting  (C.21)  for  the  hlter  h. 


(k) 

n 


'  (fc) 


M 


+  E 

il=i+l 


H 


,'(*) 

*1 


M 


+  E 


(fc)  \  hW 

~i\ii  /  *2 


*2=11+1 


(1  <  i  <  M  -  1). 

(C.24) 
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Separating  the  expanded  term  yields 

M  M  M 

*-*••  +  E  (-<’K<*,+  E  E  (-«£')  h<f  (i < * < m — i). 

*i=*+l  i%=i+ 1  *2=*i+l 

(C.25) 

Substituting  (C.21)  into  the  third  term  leads  to 

M 

hf’  =  h'<‘>+  E  (-<’) 

*1=*+1 

M  M  /  M  \ 

+  E  E  h'f>+  V  (-H«)h«  (1  <  i  <  M  —  1) 

*1=*+1  12=*1+1  \  *3=*2  +  l  / 

(C.26) 

Separating  the  last  term  on  the  right-hand  side  of  the  equation  yields 

M  M  M 

hf>  =  \{k)  +  E  (-<’K‘)  +  E  E  (-<>)  hf 

*l=*+l  *i=*+l  i2=*l+l 

M  M  M 

+  E  EE  (-<’)  0H‘H)  (-Hft)  h!f  (1  <  i  <  M  -  1)  (C.27) 

ii=i+l  *2=*i+l  *3=*2+l 

If  this  recursion  is  continued,  the  final  result  becomes 

M  M  M 

hi*>  =  h;«*>+  e  (-<) +  E  E  (-<>)(-<) 

*i=*+l  *i=i+l  *2— *i +1 

MM  M  M 

+  y  y  ...  y  y  y h^V-h^ V . . .  W.(fc) 

/  >  /  >  /  >  /  >  Y  Ul  J  \  *1*2  J  \  *M-2*M-1  J  Y  «M-1««  y  IM 

*1=*+1  *2=*1+1  *M-l=*M-2  +  l  *M=*M-1+1 

(1  <  i  <  M  -  1).  (C.28) 

The  result  of  (C.28)  can  be  used  to  also  develop  the  recursive  form  of  the  error 
variance. 


B.  EXPRESSION  FOR  THE  ERROR  VARIANCE 

Beginning  with  the  basic  definition  of  the  error  variance 

M 

^  =  -Mo)-E  hEhfb  (c.29) 
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the  summation  can  be  separated  into 


M—l 

4 = mo)  -  y  if.  (&»» 

i= 1 

Substituting  (C.28)  into  (C.30)  for  the  filter  term  representing  the  filters  for  signals 
1  through  M  —  l  leads  to 


M—l 


al  =  Rd{ 0)  -  J2  X 

i= 1 

M  MM 

E  (-i*K"+  E  E  (-<’ 


h  1 


-H^n  h  • (fc)  + 


*1*2  /  *2 


*1=*+1 


*i=*+l  *2=*i+l 
MM  M  M 

E  E  ■  E  .  E 

*1=*+1  *2=*1+1  *M-l=*M-2  +  l  lM=*M-l+l 


-Hg> )  x 


H 


(fc) 

*1*2 


H 


(fc) 


*M-2*M-1 


H 


(fc) 


h 


'(fc) 


_  ~(fc)T,  (fc)* 

ldM  11 M  ' 


(C.31) 
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The  multiple  summation  term  can  be  separated  into  two  parts  leading  to 


M—l 


=  Rd{ o)  -  i 


{k)T 

di 


X 


i=  1 


M—l 


*l=*+l 


M—l  M—l 


E  E  (-Hg?)( 

ii=i+l  *2=*i+l 


+  •  •• 


M—l  M—l  M—l  M—l 

+  E  E  E  E  (-<>)* 

*1=2+1  *2=*1+1  *M  — 1=*M— 2  + 1  *M=*M-1+1 


+  E  E  •••  E  (-<+ 

2l=2+l  22— ^1+1  —  1— 2M  —  2  +  1 


-H 


(fc) 

'*1*2 


H 


(fc) 


*M-2*M-1 


-]+^  1  h 

n>M-lM  1  n*M 


-(fc)T,  (fc)* 

X1M  • 


(C.32) 


The  first  two  terms  on  the  right-hand  side  is  the  error  variance  for  an  (M  —  1)- 
signal  system.  Thus  defining  cr|  M  =  a\  and  representing  the  error  variance  for  the 
(M  —  l)-signal  system  as  +  leads  to 


a 


2 

k,M 


=  a 


2 

k,M- 1 


( k )* 
M 


M—l 


i= 1 


M—l  M 

-e  e 

2=1  21=2+1 

M—l  MM  M 

EE  E  E  fST(-H+x 

2—  1  21=2+1  22—21+1  2jvT  — 1— ^M  — 2  +  1 


H 


( k )* 

—2'i'M  —  1 


1  h 

niif-iM  j  n*M 


'  (fc)* 


(C.33) 


143 


C.  SUMMARY 


For  clarity  and  ease  of  reference,  the  results  for  the  innovations  form  of  the 
optimal  filter  are  summarized  here.  Given  an  M-signal  multirate  system,  the  solution 
to  the  Wiener-Hopf  equations  is  given  by 


T?  (fc) 

rvn 

p  (fe) 

-1 

1 

<— 1 

•  • 

IFh 

1 

cr 

^2 

p(fc)*T 

rvliW 

p  W 

f(fc) 

1  dM 

The  following  terms  are  defined  for  use  in  the  filter  equations 


B(fc)  = 
v 


1 

>pr 

b(fc)  - 

1 - 

•• 

i 

■oj 

_ 1 

°di  ~ 

i 

and 


and 


Gik)  =  ( 


p  (k) 

'*11 

p  (k) 
**-12 

p  (k)*T 

rv12 

R?} 

-r>(k)*T  p(fc)*T 
^Ip-l)  X^2(i— 1) 


E«‘>  = 


-D  (fe) 

Kn 


i  —  1 


R 


r: 


(fc) 

'l(i-l) 

(fc) 

2(i-l) 


R 


(fc) 


-l 


Bifc)  1 <  i  <  M 


i  =  1 


R,lfc)  -  G{k)*TB{k)  1  <  i  <  M 


The  equations  for  the  filter  coefficient  in  the  innovations  representation  of  the  filter 
are  then  summarized  in  Table  C.l. 
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Table  C.l.  Innovations  Form  of  the  Optimal  Filter 
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APPENDIX  D.  MULTICHANNEL  LEVINSON 

RECURSION 


This  appendix  briefly  discusses  the  procedure  known  as  the  multichannel 
Levinson  algorithm,  (also  called  the  Levinson- Wiggins- Robinson  (LWR)  algorithm). 
For  more  detailed  discussions  of  this  topic  see  [Ref.  64,  70,  71].  The  multichannel 
Levinson  algorithm  is  a  method  of  solving  the  Normal  equations  for  multichannel  sys¬ 
tems  in  a  recursive  fashion.  This  algorithm  is  an  extension  of  the  Levinson  recursion 
for  single-channel  systems. 

A.  THE  MULTICHANNEL  SYSTEM 

In  order  to  properly  describe  the  multichannel  Levinson  algorithm,  it  is  nec¬ 
essary  to  first  consider  a  multichannel  random  process  and  its  correlation  function. 
Let  the  random  process  x[n]  have  the  form 

xi[n] 
x2[n 

xc[n]  J 

where  C  is  the  number  of  channels.  The  correlation  function  associated  with  x [n]  is 
given  by 

Rx[l]  =  £{x[n]x*T[n  —  U- 

where  R x[l]  is  a  C  x  C  matrix. 

The  multichannel  Levinson  algorithm  solves  for  the  filter  coefficients  and  pre¬ 
diction  error  covariance  matrices  for  linear  prediction  of  the  multichannel  random 


x[n  = 
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process.  In  so  doing,  it  finds  the  parameters  for  both  the  forward  and  backward 
prediction  of  the  multichannel  random  process. 

B.  GENERAL  FORM  OF  THE  ALGORITHM 

To  initialize  the  recursion,  the  following  variables  are  defined, 

R-o  =  Rx[l]  A0  =  B0  =  IcxC  S0  =  Eg  =  Rx[0], 

where  IgxC  is  an  identity  matrix  of  size  C  x  C .  In  addition,  a  cumulative  correlation 
matrix  is  defined 

Rxb] 

Rxb  - !] 

Rx[l] 

where  p  is  between  1  and  the  desired  order.  Throughout  this  algorithm,  the  variables 
with  the  superscript  b  refer  to  the  backward  prediction  parameters  and  the  other 
variables  refer  to  the  forward  prediction  parameters. 

Step  1.  Calculate  the  reflection  coefficients,  Tp  and  Tbp 

r  -|  rr  i  iT 

AP=[Rf  R^  •••  R^J  K-i  =  [[Ri  R2  •••  RPJ  b;_! 

Tp  =  p 

and 

F6  —  T_1  At 

Step  2.  Calculate  the  pth  order  prediction  coefficients,  and  Bp 
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Bp_i  0 

b;=  _  -  _ k 

0  A-p-i 

Step  3.  Calculate  the  pth  order  Error  Variance  coefficients,  and 


These  steps  are  repeated  until  the  desired  order  has  been  reached.  At  that 
point,  the  prediction  coefficients,  error  variances  and  reflection  coefficients  for  all 
orders  np  through  the  desired  order  have  been  calculated. 
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APPENDIX  E.  A  GENERALIZED  MATRIX 
TRIANGULAR  FACTORIZATION 
ALGORITHM 


In  order  to  solve  (4.41)  for  the  transform  matrix  Ao  when  developing  a  so¬ 
lution  to  the  multirate  linear  prediction  problem  in  Chapter  IV,  a  generalized  UL 
factorization  was  developed.  The  generalized  UL  factorization  presented  in  this  ap¬ 
pendix  is  modeled  after  the  standard  LU  factorization  algorithm  that  can  be  found  in 
many  linear  algebra  books,  such  as  [Ref.  65] ;  however,  when  factored,  the  left  matrix 
is  upper  triangular  and  the  right  matrix  is  lower  triangular.  The  algorithm  is  called 
“generalized”  since  it  allows  for  a  block  UL  factorization  using  different  size  blocks. 
It  thus  provides  for  Gaussian  elimination  of  multiple  rows  at  the  same  time.  This  is 
important  to  the  multirate  linear  prediction  solution,  since  eliminating  multiple  rows 
at  the  same  time  is  necessary  when  dealing  with  joint  observations  of  multiple  signals. 

A.  LU  TRIANGULAR  FACTORIZATION  REVIEW 

This  section  briefly  reviews  the  algorithm  for  computing  the  basic  LU  trian¬ 
gular  factorization  of  a  matrix.  For  a  more  detailed  discussion  of  this  methodology, 
refer  to  any  elementary  linear  algebra  textbook,  e.g,  [Ref.  65]. 

Given  an  m  x  n  matrix  A, 

flip  Ol,2  •  •  •  fll,n 

fl'2,1  «2,2  ‘  ‘  ‘  0,2, n 

flm,l  flj, 2  ‘  "  "  fl)n,n 

it  is  desired  to  decompose  the  matrix  into  upper  and  lower  triangular  matrices,  such 
that 

A  =  LU, 
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where  L  is  the  lower  triangular  matrix  with  l’s  on  the  main  diagonal  and  U  is  the 
upper  triangular  matrix.  The  procedure  for  this  method  is  as  follows: 

Step  1:  For  i  —  1,  2, . . . ,  m  —  1,  eliminate  row  i  from  rows  i  +  1  to  m: 

a.  Calculate  row  multipliers  for  rows  i  +  1  to  m 

Mkti  =  a^i/a^i  ,  k  =  i  T  1,  %  T  2, . . . ,  m 

b.  Subtract  row  i  from  rows  %  +  1  to  m 

ak,i  =  ak,i  -  Mk}i  *  aij  =  0  ,  k  =  i  +  l,i  +  2, . . .  ,  m 
ak,i  =  ak,i  ~  Mkti  *  ciij  ,  k,  l  =  i  +  1,  %  +  2, . . . ,  m 

After  eliminating  i  rows,  the  matrix  has  the  form: 

Ol,l  Ol,2  Ol,3  '  '  '  al,i  Ol,i+l  •  •  •  0,1, n 

0  02,2  O 2,3  '  '  '  02, i  02,i+l  '  •  •  02, n 

0  0  03,3  ‘  ‘  '  ^3, i  03,1+1  '  '  '  03,n 

A, 

0  0  0  *  O^i  O^, 2-|-l  '  '  '  02,21 

0  0  0  ■  Oi+i,i  Oi+i,i+i  ■  ■  ■  Oi+i,n 

0  0  0  '  On J,j  Om,2  +  l  ■  '  ■  0221,21 

Repeat  step  1  until  matrix  A  is  upper  triangular.  This  is  the  matrix  U. 


0  +  1  Oi,2  01,3  •  •  •  O  i  ,2  Oi,  j+i  •  •  •  Oi,2i 

0  02,2  02,3  •  •  •  02,2  02,2+1  '  '  '  02,21 

0  0  03,3  •  •  •  03,2  03,2-1-1  •  •  •  03,21 

U  =  A„i  = 

0  0  0  *  0^,2  02,2+1  *  *  *  02,21 

0  0  0  ■  0  Oi-f-1 ,2  +  1  '  '  '  02-1-1,21 

0  0  0  ■  ■  •  0  0  •  •  •  a,2i,n 
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Step  2:  The  lower  triangular  matrix  can  be  formed  as 


L 


1  0 

M2,  i  1 

M:iA  M3j  2 


0 

0 

1 


0 

0 

0  , 


Mm>  3  •  •  •  1 

where  Mt]  are  defined  above. 


B.  A  GENERALIZED  UL  TRIANGULAR  FACTORIZA¬ 
TION  ALGORITHM 

The  basic  LU  triangular  factorization  algorithm  can  be  modified  and  extended 
to  a  generalized  UL  triangular  factorization  algorithm.  This  new  procedure  requires 
the  defining  and  use  of  interim  variables  to  transform  an  original  matrix  A. 

Given  an  m  x  n  matrix  A. 

Ol,2  •  •  •  « 1  ,n 

0*2,2  ‘  ‘  «2 ,n 

2  *  *  *  ^m,n 

it  is  desired  to  decompose  the  matrix  into  generalized  upper  and  lower  triangular 
matrices,  such  that 

A  =  UL, 

where  U  is  a  block  upper  triangular  matrix  and  L  is  a  block  lower  triangular  matrix. 
The  blocks  are  not  all  of  the  same  size,  but  the  sizes  of  the  desired  blocks  are  known 
(or  given). 

The  procedure  for  this  method  is  as  follows  (starting  with  the  last  column ), 


A  = 


«i,i 

«2,1 

1 
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Step  1:  Eliminate  columns  i  —  r  +  1  through  i  from  columns  1  through  i  —  r,  where  r 
is  the  number  of  columns  to  be  eliminated  at  the  same  time  and  i  is  the  index  of  the 
last  of  these  columns. 

a.  Define  the  following  vectors  and  matrices 


&i— r+1, 2— r+1 


A 


■(i-r+l)t,(i-r+l)i 


r+1 


r+l,z 


^ i,i 


a(i-r+l)i,k 


,  k  =  1, 2, . . . , i  —  r 


®i,k 


al,(i—r+l)i 


i—r+1 


1,2, 


!  1 


—  r 


b.  Calculate  column  multipliers  for  columns  1  to  i  —  r 

M(i-r+\)i,k  =  ^-(i_r+l)i,(i_r+1)ia(i-r+l)i,fci  ^  =  1,  2,  .  .  .  ,  i  —  V 

c.  Subtract  columns  i  —  r  +  1  through  i  from  columnss  1  through  i  —  r 


a(i— r+l)i,fc  a(i— r+l)i,fc  -A-(i — — r-+l)i,fc  0  ,  A  1,  2,  .  .  .  ,  i  T 

&l,k  +,fc  aZ,(i— r+l)jAf(j_ r+l)i,fc  j  k,l  1,  2,  .  .  .  ,  i  f 


At  the  ith  iteration,  the  matrix  has  the  form: 
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fll,l 

^1  ,i—r 

. .  T 
- 

+ 

Al,n— 2 

&l,n—  1 

&l,n 

a3,l 

Qj i — r,i — r 

ai-r,(i-r+l)i 

flj— r,n—  2 

flj— r+l,n— 2 

r,n— 1 

r+l,n—  1 

^2— r,n 

r+l,n 

aip— r 

a(i-r+l)i,i-r 

Ap_ r+l)i,(i— r+l)i 

flj,n— 2 

^i,n 

0 

0 

0 

0  •• 

Am— 2,m— 2 

&m—2,n—l 

Q"m— 2,n 

0 

0 

0 

0  •• 

0 

Q"m— l,n— 1 

(^m— l,n 

0 

0 

0 

0  •• 

0 

0 

Qjm,n 

Repeat  step  1  until  matrix  A  is  transformed  into  the  upper  triangular  matrix  U. 


flip 

fll,2 

fll,3  •  • 

flip 

fllp+1 

^l,n 

0 

fl'2,2 

fl2,3  ‘  ‘ 

fl'2p 

fl2p+l 

^2,n 

0 

0 

fl3,3  •  • 

A  3.7 

fl3p+l 

^3,n 

U  =  Am  = 

0 

0 

0  •• 

flip 

Ajp+l 

^i,n 

0 

0 

0  •• 

•  0 

Ai+lp+1 

0 

0 

0  •• 

•  0 

0 

Q"m,n 

Step  2:  After  the  2nd  row  is  eliminated  from  the  1st  row,  the  upper  and  lower  trian¬ 
gular  matrices  can  be  formed  as 


L 


1  0  0 

M2,i  1  0 

M3ii  M3)2  1 


0 

0 

0 


MTOji  Mto,2  MTOj3  •  •  •  1 
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C.  EXAMPLE 


Given  a  matrix  A, 


A  = 


fll,l 

:  0  i  .2 

Ol,3  : 

:  <2-1, 4  : 

:  01,5 

Ol,6 

<3-2,1 

:  ®2,2 

0-2,3  : 

:  02,4  : 

:  °2,5 

0-2,6 

03,1 

:  03,2 

03,3  : 

:  <23,4  : 

:  03,5 

03,6 

04,1 

04,2 

04,3  : 

!  4  ! 

:  04,5 

04,6 

<3-5,1 

0-5,2 

05,3 

05,4  : 

:  05,5 

05,6 

06,1 

06,2 

06,3 

06,4  : 

:  Og,5 

^6,6 

(E.l) 


where  at^  is  the  matrix  element  for  the  ith  row  and  jth  column.  It  is  desired  to  decom¬ 
pose  the  matrix  using  the  generalized  UL  algorithm.  For  this  particular  example,  it 
also  is  desired  to  eliminate  columns  5  and  6  together  and  followed  by  the  elimination 
of  columns  2  and  3  together,  as  shown  by  the  partitioning  in  (E.l).  The  number  of 


12  12 


,  where 


columns  to  be  eliminated  can  be  represented  in  vector  form  as 
the  numbers  refer  to  the  block  sizes  (number  of  columns  to  be  eliminated  together). 
Step  1:  Eliminate  columns  5  and  6 

a.  Define  the  following  vectors  and  matrices 
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As6,56 


05,5  ®5,6 

06,5  «6,6 


05,1 

05,2 

05,3 

05,4 

a56,l  — 

a56,2  — 

a56,3  — 

a56,4  — 

06,1 

06,2 

06,3 

06,4 

al,56  — 


Ol,5  al,6 


a2,56  — 


®2,5  02,6 


a3,56  — 


®3,5  03,6 


a4,56  — 


04,5  04,6 


b.  The  multipliers  for  columns  1  through  4  are 


M56ji  —  Agg'^ga.sfj,  1  M56)2  —  Agg'^ga.sfj^ 

M56,3  =  Agg'gga.sfj^s  M56)4  =  Agg'gga.sG^ 

c.  Eliminating  columns  5  and  6  from  columns  1  through  4  produces  (primes 
have  been  added  to  the  modified  values,  for  this  example  only,  to  highlight  the  changes 
at  each  step) 


r 

■56,1  = 

a56,l  — 

A56, 56^56,1  = 

°1,1 

=  Oi,i 

—  al, 56^56,1 

°3,1 

=  03,1 

—  a3, 56^56,1 

/ 

■56,2  = 

a56,2  — 

A56,56M56,2  = 

r 

<3-1,2 

=  di,2 

—  al, 56^56, 2 

r 

a3,2 

=  03,2 

—  a3, 56^56, 2 

/ 

-56,3  = 

a56,3  — 

As6,56  AE5g,3  = 

al,3 

=  Oi,3 

—  al,  56^56, 3 

r 

a3,3 

=  03,3 

—  a3,56  AE5g,3 

/ 

■56,4  = 

a56,4  ~ 

A56,56M56,4  = 

r 

al,4 

=  ai,4 

~  al, 56^56, 4 

r 

a3,4 

=  03,4 

—  a3, 56^56, 4 

a2,l  —  <3-2,1  —  a2,56  AE561 

«4,1  =  04,1  —  a4, 56^156,1 

°2,2  =  fl2,2  —  a2, 56^156,2 
0-4, 1  =  04,2  ~  a4, 56^156,2 

a2,3  =  02,3  —  a2,56  AE5g,3 

O43  =  04,3  —  a4,5gM56,3 

0-2,4  =  0-2,4  —  a2,56  ^[55,4 

04,4  =  04,4  —  a4,56  AE5g,4 
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d.  Matrix  A  with  columns  5  and  6  eliminated  looks  like 


al,l 

.  / 

:  a12 

al,3 

.  / 

:  0 1  4 

:  al,5 

Ol,6 

1 

a2,l 

•  / 

:  a22 

1 

°2,3 

•  / 

•  °2,4 

:  a2,5 

a  2,6 

a3,l 

.  / 

:  a3,2 

1 

a3,3 

.  / 

:  °3,4 

:  a3,5 

03,6 

1 

a4,l 

/ 

a4,2 

/ 

a4,3 

•  / 

(3-4  4 

:  «4,5 

04,6 

0 

0 

0 

0 

:  a5,5 

05,6 

0 

0 

0 

0 

:  a6,5 

«6,6 

Step  2:  Eliminate  row  4 

a.  The  multipliers  for  columns  1  through  3  are 

M4,i  =  Cl^\/ U44  M4i2  —  Cl4  2 / ®4)4  M4i3  =  04  3/ a44 


b.  Eliminating  column  4  from  columns  1  through  3  produces 


11 

°4,1  = 

/ 

a4,l  — 

®4i4d'h4,l  — 

II 

a2,l 

=  a2,l 

®2,4'^/^4,l 

II 

a4,2  = 

/ 

a4,2  — 

®4i4T/4i2  = 

II 

a2,2 

r 

=  a2,2 

Q'2,4-^-^4,2 

II 

°4,3  = 

/ 

°4,3  — 

a44M4)3  = 

II 

°2,3 

r 

=  a2,3 

®2,4^4,3 

//  /  /  ,  , 

&3  i  —  (X3  ^  (3^  4^^4,1 

//  /  /  ,  , 

1  —  o-i  ^  (3i  4-^41 

//  /  '  T\  /T 

a3,2  =  a3,2  —  0t3,4-'“4,2 
//  /  /  71  /T 

2  ~  ^1  2  ^1  4-^*4, 2 

//  /  '  n  /r 

a3,3  =  a3,3  —  a3,4-^4,3 

//  /  I  ,  /r 

(2i3  =  ^1,3  ^'1,4'^4,3 
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c.  Matrix  A  with  columns  4  through  6  eliminated  looks  like 


//.////./ 


'1,1  ■ 

:  al2 

al,3  ' 

:  a14  : 

:  °1,5 

d  1,6 

// 

■  n 

// 

■  r 

2,1  ' 

■  a2,2 

a2,3  ' 

:  o24  : 

:  a2,5 

02,6 

// 

■  U 

H 

■  f 

'3,1  : 

'■  °3,2 

a3,3  : 

:  a3,4  : 

:  a3,5 

03,6 

0 

0 

0  : 

/ 

CL  4  4 

:  «4,5 

04,6 

0 

0 

0 

0  : 

:  a5,5 

05,6 

0 

0 

0 

0 

'■  a6,5 

06,6 

Step  3:  Eliminate  columns  2  and  3 

a.  Define  the  following  vectors  and  matrices 


^23,23  — 

// 

a2,2 

// 

// 

a2,3 

// 

3-23,1  — 

// 

a2,l 

n 

al,23  ~ 

//  // 
°1,2  al,3 

a3,2 

a3,3 

a3,l 

b.  The  multiplier  for  column  1  is 

M23,i  =  A231j23a23)i 

c.  Eliminating  columns  2  and  3  from  column  1  produces 


*23,1 


—  ^23,1  —  A23  93M23  1  —  0 


‘1,1 


=  a 


1,1 


AL,23M 


23,1 
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d.  Matrix  A  with  columns  2  and  3  eliminated  looks  like 


III 

°1,1 

.  n 

:  a12 

II 

al,3 

.  / 

:  a  j  4 

:  °1,5 

a  i,6 

0 

•  II 

•  a2,2 

II 

a2,3 

•  / 

•  °2,4 

:  a2,5 

0'2,6 

0 

•  II 

’■  °3,2 

II 

a3,3 

.  / 

’  °3,4 

:  a3,5 

^3,6 

0 

0 

0 

•  / 
a  4  4 

:  «4,5 

04,6 

0 

0 

0 

0 

:  a5,5 

05,6 

0 

0 

0 

0 

'■  a6,5 

^6,6 

Step  4:  Build  the  U  and  L  matrices 

a.  The  matrix  U  is  the  transformed  matrix  A 


/// 

°i,i 

.  // 

:  a1 2 

II 

al,3 

.  / 

:  a14 

:  al,5 

Ol,6 

0 

•  // 

•  a2,2 

II 

0'2,3 

•  / 

•  ®2,4 

'■  a2,5 

02,6 

0 

•  // 

:  a3,2 

II 

a3,3 

.  / 

:  a3,4 

'■  03j5 

03,6 

0 

0 

0 

•  / 

&4  4 

:  O45 

04,6 

0 

0 

0 

0 

:  05j5 

05,6 

0 

0 

0 

0 

:  06j5 

^6,6 
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b.  The  matrix  L  is  made  up  of  the  multipliers,  such  that 

r  i  :  o  o  ;  o  ;  o  o 

i  o  ;  o  ;  o  o 

o  i  i  o  ;  o  o 

d'^4.2  Ma,3  •  1  •  0  0 

•10 

M56)i  M56i2  M56)3  M56i4 

:  0  1 


M23,l 

L  = 

M41 
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